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I.  FOREWARD 


The  research  program  has  been  aimed  at  developing  computational  procedures 
for  simulation  of  flow  about  projectiles.  During  the  three  year  effort  a  variety  of 
ideas  were  pursued  and  implemented.  As  a  result  this  final  report  will  include 
summaries  on  numerical  algorithm  developments,  methods  of  grid  generation,  and 
descriptions  of  flow  field  solution  codes  for  projectiles  in  the  transonic  range.  Some 
five  technical  papers  have  been  published  which  describe  much  of  the  research  effort. 
These  are  included  in  this  final  report  and  they  constitute  the  technical  content  of 
this  report.  Some  ongoing  work  which  has  not  yet  reached  fruition  is  also  described. 

The  research  effort  was  initially  directed  toward  two  tasks:  1)  basic  algorithm 
developments;  and  2)  development  of  a  parabolized  Navier-Stokes  (PNS)  computer 
code  to  solve  finned  projectiles  in  supersonic  viscous  flow.  The  second  task  was 
later  redirected  to  the  development  of  computer  codes  for  computing  the  transonic 
flow  about  projectiles  with  base.  This  effort  was  undertaken  in  collaboration  with 
Messrs.  Nietubicz  and  Sahu  of  the  Ballistics  Research  Laboratory  (BRL).  (The 
original  task  of  computing  the  supersonic  flow  about  finned  projectiles  using  a  PNS 
code  was  subsequently  accomplished  by  Mann  Mohan  Raj,  et.  al.,  under  joint  BRL 
and  NASA  Ames  Research  Center  sponsorship.) 

The  research  effort  has  greatly  benefited  by  a  sustained  collaboration  with  re¬ 
searchers  at  BRL.  This  collaboration  properly  focused  the  work  on  realistic  prob¬ 
lems,  help  stimulate  new  concepts,  and  provided  necessary  stimulate  new  concepts, 
and  provided  necessary  critiques. 

The  report  is  divided  into  four  main  areas.  Section  II  describes  the  algorithms 
and  methodology  for  computing  transonic  flow  about  projectiles  with  base.  Section 
III  describes  some  work  in  grid  generation,  while  Section  IV  contains  a  potpourri 
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II.  TRANSONIC  PROJECTILES  WITH  BASE  FLOW 

A  major  accomplishment  of  the  research  program  has  been  the  development  of 
computer  codes  to  simulate  transonic  projectiles  with  base  flow.  This  project  was 
carried  out  in  collaboration  with  J.  Sahu  and  C.  J.  Nietubicz  at  BRL  with  Sahu 
being  responsible  for  the  code  implementation. 

A  full  description  of  the  numerical  procedures  is  given  in  the  two  appended 
ALAA  papers.  A  key  feature  of  the  projectile  base  flow  code  is  its  segmentation 
concept.  Beginning  with  a  basic  projectile  and  sting  code  due  to  Nietubicz,  a  simple 
way  of  dividing  up  the  computational  domain  was  devised  which  maintained  the 
simplicity  of  the  implicit  numerical  algorithm.  The  sketches  shown  in  Figure  2 
of  the  second  attached  ALAA  paper  (83-0224)  illustrate  the  segmentation  process. 
The  idea  here  is  that  we  add  grids  as  needed  yet  solve  the  flow  field  implicitly  as 
one  large  grid  in  computational  space.  Various  flags  are  used  to  properly  turn-off 
or  connect  the  domains  together.  This  concept  can  be  extended  to  more  complex 
geometries,  and  currently,  Sahu  is  treating  a  projectile  with  a  cut-out  in  the  base 
by  using  three  grid  segments. 
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Abstract 

The  Azimuthal- Invariant  Thin-Layer  Navier- 
Stokes  computational  technique  has  been  modified 
for  projectile  base  flow  analysis.  The  resulting 
new  numerical  capability  is  used  to  compute  the 
entire  projectile  flow  field  including  the  recir- 
culatory  base  flow.  Computed  results  show  the 
qualitative  and  quantitative  details  of  the  over¬ 
all  base  flow  structure.  Base  drag  is  computed 
for  a  secant-ogive-cylinder  projectile  and  com¬ 
pared  with  limited  experimental  and  semi-empirical 
data.  Results  are  also  presented  which  show  the 
variation  of  pressure  drag,  skin  friction  drag  and 
the  total  aerodynamic  drag  for  Mach  numbers  .9  <  M 
<  1.2. 

Nomenclature 

a  speed  of  sound 

A  cross  sectional  area  at  the  base 

C|)b  base  drag  coefficient,  2  0b/p—u£/t 

Cp  specific  heat  at  constant  pressure 

Cp  pressure  coefficient,  2(p  -  pj/p_a2 

D  body  diameter  (57.15mm) 

Db  base  drag 

e  total  energy  per  unit  volume/p^a2 

q  vector  of  dependant  variables  in 

.  .  transformed  equations 

E,  F  flux  vector  of  transformed  Navier-Stokes 
.  equations 

H  n-invariant  source  vector 

J  Jacobian  of  transformation 

M  Mach  number 

p  pressure/p„a2 

Pr  Prandtl  number, 

R  body  radius  p 

Re  Reynolds  number, 

S  viscous  flux  vector 

t  physical  time 

u,v,w  Cartesian  velocity  components/a,. 

U,V,W  Contravariant  velocity  components/a,,. 
x,y,z  physical  Cartesian  coordinates 
a  angle  of  attack 

r  ratio  of  specific  heats 

x  coefficient  of  thermal  conductivity 

u  coefficient  of  viscosity 

t,n,c  transformed  coordinates  in  axial, 

circumferential  and  radial  directions 
p  density/ p„ 

t  transformed  time 
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0  circumferential  angle 

Superscript 

*  critical  value 

Subscript 

b  base 

p  pressure 

v  viscous 

“  free  stream  conditions 


I.  Introduction 

The  rising  costs  of  experimental  measurements 
has  resulted  in  alternate  means  of  determining  the 
aerodynamics  of  shells.  Because  of  the  recent 
advances  in  computer  processors,  numerical  compu¬ 
tational  capabilities  have  been  developed  to  pre¬ 
dict  the  aerodynamic  behavior  of  artillery  shells. 

Recent  papers1’2  have  reported  the  develop¬ 
ment  and  application  of  the  Azimuthal-Invariant 
Thin-Layer  Navier-Stokes  computational  technique 
to  predict  the  flow  about  slender  bodies  of  revo¬ 
lution  at  transonic  speeds.  References  1  and  2 
showed  the  technique  to  be  a  viable  computational 
tool  for  predicting  both  external  and  internal 
flows  for  spinning  and  non-spinning  bodies  of 
various  geometric  shapes.  The  base  flow  of  the 
projectiles  however,  was  not  computed.  Instead 
the  projectile  base  was  modeled  by  an  extended 
sting.  Experimental  base  flow  data  is  difficult 
to  obtain  and  therefore  only  limited  data  is 
available.  No  sophisticated  numerical  techniques 
have  yet  been  utilized  for  base  flow  of  projec¬ 
tiles  at  transonic  speeds.  The  objective  of  this 
research  was  to  develop  a  new  numerical  capability 
to  compute  the  flow  field  in  the  base  region  of 
projectiles  at  transonic  speeds  and  be  able  to 
compute  the  total  aerodynamic  drag. 

The  total  drag  for  projectiles  can  he  divided 
into  three  components:  (1)  pressure  (wave)  drag; 

(2)  viscous  (skin  friction)  drag;  and  (3)  base 
drag.  For  a  typical  shell  at  H  =  .9,  the  relative 
magnitudes  of  the  aerodynamic  drag  components  are: 
(I)  pressure  drag,  20%;  (2)  viscous  drag,  30%;  and 

(3)  base  drag,  50%.  In  order  to  predict  the  total 
drag  for  projectiles,  computation  of  the  full  flow 
field  (including  the  base  flow)  must  be  made. 
Computation  of  base  flow  is  especially  important 
at  transonic  speeds. 


The  critical  aerodynamic  behavior  of  projec¬ 
tiles  occurs  in  the  transonic  speed  regime.  This 
can  be  attributed  to  the  complex  shock  structure 
which  exists  for  the  projectiles  at  transonic 
speeds.  Figure  1  is  a  spark  shadowgraph  which 
shows  the  shock  structure  for  a  typical  projectile 
at  M  .9$,  a  -  I).  It  also  shows  a  clearly 
defined  wake  behind  the  base  of  the  projectile 
devoid  of  any  vortex  shedding.  Primary  emphasis 
is  focused  on  the  base  region  flow  field  computa¬ 
tions;  however,  the  technique  used  computes  the 
full  flow  field  over  the  projectile  (including  the 
base  region).  Therefore,  all  three  components  of 
the  drag  are  computed. 

A  brief  description  of  the  governing  equa¬ 
tions  and  the  method  of  solution  are  given  in 
Sections  II  and  III.  A  unique  flow  field  segmen¬ 
tation  procedure  and  the  implementation  of  bound¬ 
ary  conditions  are  discussed  in  Section  III.  In 
Section  IV  computed  results  are  given  for  tran¬ 
sonic  flow  about  a  S-caliber  secant-ogi ve-cy) inder 
shape  for  .9  <  M  <  1.2,  u  =  0.  Velocity  vector 
plots  and  stream  function  contour  plots  are  pre¬ 
sented  to  show  the  qualitative  features  of  the 
flow  field  in  the  base  region.  AH  three  compon¬ 
ents  of  drag  are  obtained.  Base  drag  is  compared 
with  experimental  and  semi -empirical  data  while 
the  total  drag  is  compared  with  the  only  available 
semi-empirical  data.  The  encouraging  results  show 
that  the  present  computational  technique  can  be 
successfully  used  to  predict  the  base  region  flow 
field  of  projectiles.  Although  results  here  are 
reported  for  transonic  speeds,  future  computation¬ 
al  efforts  will  be  directed  at  supersonic 
velocities. 


II.  Governing  Equations 

The  Azimuthal  Invariant  (or  Generalized  Axi- 
symmetric)  thin-layer  Navier-Stokes  equations  for 
general  spatial  coordinates  4,  n,  c  can  be  written 
as1 

a  q  *  if  +  a  G  ♦  H  =  Re'* 8S  (1) 

*  k*  S 


pW 

pvw+cyp 
G  -  J~'  pvW+CyP 
»wW+ C^o 
(etp)W-r,tp 


0 

0 

H  J'1*,  uV[R.(ll-f,t)  *  «,(W-,t)] 
-pVR*n(V-nt)  -  p/(R*n) 

0 


The  thin  layer  viscous  terms  valid  for  high 
Reynolds  number  flow  are  contained  in  the  vector 

S,  where 


S  -- 


m(  WS)UC  *  (u/3)<  Vr/;yVVcKx 
“( VV  4z)wc  +  (u/3)<  WcyYczw;);z 

,  2  2  2  2  2  2  i 

t(yyy[°-^(u  +v  +  <Pr 

(  y-1  )'*  (a2)^]  +  (u/3)(cxu+cyv-*czw) 

'vvvvyyf 


The  velocities 


where  4  *  4(x,y,z,t)  is  the  longitudinal 
coordinate 

n  =  n(y,z,t)  is  the  circumferential 
coordinate 


d  -  4t  +  4xu  +  4 yv  +  4?w 
V  =  nt  +  nxu  +  nyv  +  nzw  (2) 


4  *  c(x,y,z,t)  is  the  near  normal 
coordinate 

i  =  t  is  the  time 

The  notation  for  the  physical  coordinates  x,  y,  z, 
and  the  transformed  coordinates  5,  n,  <  are  shown 

in  Figure  2.  The  vector  of  dependent  variables  q 


represent  the  contravariant  velocity  components. 

The  Cartesian  velocity  components  (u,v,w)  are 
non dimensional ized  with  respect  to  (the  free 
stream  speed  of  sound).  The  density  (p)  is  refer¬ 
enced  to  p^  and  total  energy  (e)  to  p„a£.  The 
local  pressure  is  determined  using  the  equation  of 
state. 


and  the  flux  vectors  E,  G,  H  are  given  as 


P  -  (Y-l)[e  -  0.5p(u2»v2+w2)]  (3) 


p 

pU 

PU 

puU*  Cxp 

pv 

,  E  - 

evil*  f,yp 

PW 

ewU*  '.zp 

e 

(e^p)iJ- ;.tp 

where  y  is  the  ratio  of  specific  heats. 

In  high  Reynolds  number  flows  the  thin-layer 
approximation  is  often  used  because,  due  to  com¬ 
puter  speed  and  storage  limitations,  fine  grid 
spacing  can  only  be  provided  in  one  coordinate 
direction.  The  grid  spacing  available  in  other 
directions  ts  usually  too  coarse  to  resolve  the 
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viscous  terms.  Essentially,  all  the  viscous  terms 
in  the  coordinate  direction  C  and  n  are  neglected 
while  terms  in  the  near  normal  direction  to  the 
body  c  are  retained.  The  thin-layer  generalized 
axisymmetric  equations  (1)  are  obtained  from  the 
three  dimensional  equations  by  making  use  of  two 
restrictions:  (1)  all  body  geometries  are  of  an 

axi  symmetr ic  type;  and  (ii)  the  state  variables 
and  the  contravari ant  velocities  do  not  vary  in 
the  circumferential  direction  (n).  Essentially, 
the  n-der i vat i ve  term  in  the  three  dimensional 

equations  is  replaced  by  a  source  term  H  as  it 
appears  in  equation  (1).  The  details  can  be  found 
in  Reference  1  and  ?. 

Equation  (1)  contains  only  two  spatial  deriv¬ 
atives;  however  it  retains  all  three  momentum 
equations  thus  allowing  a  degree  of  generality 
over  the  standard  axisymmetric  equations.  In 
particular,  the  circumferential  velocity  is  not 
assumed  to  be  zero  thus  allowing  computations  for 
spinning  projectiles  or  swirl  flow  to  be 
accomplished. 

For  the  computation  of  turbulent  flows  a  tur¬ 
bulence  model  must  be  supplied.  In  the  present 
calculations  a  Cebeci-type  two  layer  algebraic 
eddy  viscosity  model  as  modified  by  Baldwin  and 
Lomax3  is  used.  In  their  two  layer  model  the 
inner  region  follows  the  Prandtl-Van  Oriest  formu¬ 
lation.  Their  outer  formulation  can  be  used  in 
wakes  as  well  as  in  attached  and  separated  bound¬ 
ary  layers.  In  both  the  inner  and  outer  formula¬ 
tions  the  distribution  of  vorticity  is  used  to 
determine  length  scales  thereby  avoiding  the 
necessity  of  finding  the  outer  edge  of  the  bound¬ 
ary  layer  (or  wake).  The  magnitude  of  the  local 
vorticity  for  the  axisymmetric  formulation  is 
given  by 

It  should  be  noted  that  the  turbulence  model 
has  not  been  tailored  for  use  in  base  flow 
regions.  Moreover,  the  no  slip  boundary  condition 
is  not  applied  at  the  projectile  base  and  slip  is 
allowed  along  the  base  (inviscid  boundary  condi¬ 
tion).  The  velocity  component  normal  to  the  base 
is  however  set  to  zero. 


III.  Nunerical  Method 

a.  Computational  Algorithm 

An  implicit  approximate  factorization 
finite-difference  scheme  in  delta  form  is  used  as 
described  by  Beam  and  Warming1*.  An  implicit 
method  was  chosen  because  it  permits  a  time  step 
much  greater  than  that  allowed  by  explicit 
schemes.  For  problems  in  which  the  transient 
solution  is  of  no  interest,  this  offers  the  pos¬ 
sible  advantage  of  being  able  to  reach  the  steady 
state  solution  faster  than  existing  explicit 
schemes. 


reduces  the  solution  process  tu  one-dimensional 
problems  at  a  given  time  level.  Central  differ¬ 
ence  operators  are  employed  and  the  algorithm  pro¬ 
duces  block  tridiagonal  systems  for  each  space 
coordinate.  The  main  comput a t  i ona  1  work  is  con¬ 
tained  in  the  solution  of  these  block  tridiagonal 
systems  of  equations. 

b.  f  l  nite  Difference  i  guat  ions. 

The  resultinq  finite  difference 
equations,  written  in  delta  form  are 

( 1*  h*.An-t  jj~  *V, AtJ) ( I*  hf  ,Cn-c  J'  •  v-  t.  .1 

•  hRe"  Jo,.r  *MnJ)  *  (qn+1-qn)  -At(  a.En*  v  .fi"  (b) 

-Re'16in)-atHn-cEJ-i[(v^f)z  ,  [v^kjJqH 


Here  h  =  ut  because  only  first  order  accuracy  in 
the  time  differencing  is  needed  for  the  steady 
state  flows  which  are  considered  here.  This 
choice  corresponds  to  the  Euler  implicit  time 

differencing.  The  6’s  represent  central  differ¬ 
ence  operators,  a  and  V  are  forward  and  backward 
difference  operators  respectively.  The  Jacobian 

*  dE  <$G 

matrices  A  =  — ,  C  =  — =  along  with  the  coeffi- 

aq  3q 

cient  matrix  M  obtained  from  the  local  time 

linearization  of  S  are  described  in  detail  in 
Reference  6.  Fourth  order  explicit  (t^)  and 

implicit  ( e j )  numerical  dissipation  terms  are 

incorporated  into  the  differencing  scheme  to  damp 
high  frequency  growth  and  thus  to  control  the 

nonlinear  instabilities.  A  typical  range  for  the 
smoothing  coefficients  is  e-  =  (1  to  5)  At  with 
£j  «  3eE. 

c.  Flow  Field  Segmentation 

The  objective  is  to  compute  the  full  flow 
field  (including  the  base  region)  of  a  projectile 
at  transonic  speeds.  Figure  3  shows  a  schematic 
illustration  of  the  flow  field  segmentation  used 
in  this  study  for  computational  purposes.  The 
hatched  region  represents  the  projectile.  The 
region  ABCD  includes  the  projectile  base  and  the 
wake  and  will  be  referred  to  as  the  base  region. 

The  curvilinear  coordinates  used  for  the 
longitudinal  and  normal  directions  are  represented 
by  their  indices  0  and  L.  The  line  J  =  1  starts 
at  the  downstream  boundary  (line  CO)  in  the  base 
region.  J  is  incremented  until  the  line  J  *  JB  is 
reached  which  represents  the  base  of  the  projec¬ 
tile.  The  line  J  *  JB+1  is  at  the  nose  of  the 
projectile  and  J  is  then  incremented  until  the 
line  0  *  JMAX  is  reached  which  is  the  downstream 
boundary  in  the  outer  region. 


The  Beam-Warming  impl ici t_  algori thm  has 
been  used  in  various  appl  ications l" 6.  The  algo¬ 
rithm  can  be  first  or  second  order  accurate  in 
time  and  second  or  fourth  order  accurate  in  space. 
The  equations  are  factored  (spatially  split)  which 


As  for  the  other  coordinate  L,  in  the  base 
region  L  *  1  starts  at  line  AC  which  is  a  computa¬ 
tional  cut  through  the  physical  wake  region.  L  is 
incremented  until  L  -  LMAX  which  is  the  line  of 
symmetry  (line  BO).  In  the  outer  region  l  -  1 
starts  out  from  the  projectile  surface  and  L  is 


incremented  all  the  way  to  the  outer  boundary 
where  L  =  LMAX. 

Implicit  integration  is  carried  out  from  J  - 
2  to  J  «  JB-1  and  from  J  =  JB+t  to  J  =  JMAX-l  in 
the  longitudinal  direction  and  from  L  =  2  to  L  = 
LMAX-1  in  the  normal  direction. 


Implementation  of  Boundary  Conditions 


The  no  slip  boundary  conditions  for 
viscous  flow  is  enforced  by  setting 


•j  =  v  w  o 


(6) 


on  the  projectile  surface  except  for  the  base.  At 
the  projectile  base  (J  *  JB)  the  velocity  compon¬ 
ent  normal  to  toe  base  is  set  to  zero,  i.e.  U  *  0, 
while  other  flow  variables  are  set  to  be  equal  to 
those  at  J  =  JB-l.  In  other  words,  slip  is  allow¬ 
ed  along  the  base  (inviscid  boundary  condition). 
Future  work  will  be  directed  at  the  implementation 
of  viscous  boundary  condition  at  the  base  to 
further  access  this  approximation. 

Care  must  be  taken  in  the  implementation  of 
the  boundary  conditions  along  line  AC  which  is  the 
computational  cut.  After  trial  and  error  the  flow 
variables  above  and  below  the  cut  were  simply 
averaged  to  determine  the  boundary  conditions  on 
the  cut.  This  procedure  proved  to  work  well.  On 
the  centerline  of  the  wake  region,  a  symmetry 
condition  is  imposed. 


3v 

■J z 
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Free  stream  conditions  are  used  at  the  outer 
boundary.  Simple  extrapolation  for  all  flow  vari¬ 
ables  is  used  at  the  downstream  boundary  (lines  J 
»  1,  JMAX) .  During  transient  calculations,  use  of 
a  specified  outflow  pressure  can  give  rise  to 
numerical  oscillations  in  the  base  region  flow 
field.  Eventually,  these  grow  and  swamp  the  solu¬ 
tion.  This  difficulty  is  avoided  by  simply  extra¬ 
polating  pressure  to  the  downstream  boundary  which 
is  the  procedure  always  used  with  supersonic  out¬ 
flow.  A  combination  of  extrapolation  and  symmetry 
is  used  at  J  *  JB+ 1 . 

As  a  result  of  the  flow  field  segmentation 
procedure  described  in  Section  III  b,  the  block 
tridiagonal  matrix  in  the  t  direction  has  elements 
at  J  »  JB,  JM  which  are  treated  as  internal 
boundaries  in  the  computational  domain.  The  block 
tridiagonal  matrix  in  the  f,  direction  takes  the 
following  form  (after  setting  Cj  =  0  to  simplify 
the  illustration) 


Here  A's  denote  the  qua  . -.j  A  and  I  is  a  5x5 

identity  matrix.  No  ;he  appearance  of  two 
uncoupled  block  tridia.  ,  The  rows  at  JB  and 
JB+1  are  particularly  y  e  because  boundary 
conditions  are  updated  >.  ,-iicitly  at  the  end  of 
inversions.  These  changes  were  easily  implemented 
in  a  modular  fashion  into  an  existing  code  for 
projectile  base  flow  computations.  One  simply 
fills  the  block  tridiagonal  matrix  ignoring  the 
base  JB  and  the  nose  axis  JM.  Elements  in  these 
rows  are  then  overloaded  as  shown  above.  The  flow 
field  segmentation  does  not  affect  the  block  tri¬ 
diagonal  matrix  in  the  i  direction. 


IV.  Results 

A  series  of  computations  have  been  made  for 
the  3  caliber  (1  caliber  =  1  max.  body  diameter) 
secant-ogive  nose  and  3  caliber  cylinder  shape 
shown  in  Figure  4.  All  the  computations  were 
obtained  for  Mach  numbers  .9  <  M  <  1.2  and  a  =  0. 
Limited  experimental  base  pressure  measurements 
have  been  made  by  Kayser7  for  this  projectile 
shape  and  compared  with  the  computed  results.  The 
projectile  base  was  supported  by  a  sting  attached 
to  it  and  meaurements  of  base  pressure  were  made 
at  only  one  location  along  the  base.  These 
experiments  were  conducted  at  Langley  Research 
Center  8- foot  Transonic  Pressure  Tunnel.  Compu¬ 
tational  base  pressure  results  are  also  compared 
with  available  semi -empirical11  data.  The  results 
are  presented  in  the  form  of  surface  pressure 
distribution,  contour  plots  and  velocity  vector 
plots. 

The  computational  grid  used  for  the  numerical 
computations  was  obtained  from  a  versatile  grid 
generator  developed  in  Reference  9.  This  program 
allows  arbitrary  grid  point  clustering,  thus  ena¬ 
bling  grid  points  for  the  projectile  shapes  to  be 
clustered  in  the  vicinity  of  the  body  surface. 
The  grid  consists  of  108  points  in  the  longitudi¬ 
nal  direction  and  50  points  in  the  radial  direc¬ 
tion.  The  full  grid  is  shown  in  Figure  5  while 
Figure  6  shows  an  expanded  view  of  the  grid  in  the 
vicinity  of  he  projectile.  The  computational 
domain  extended  to  4  body  lengths  in  front,  4  body 
lengths  in  the  radial  direction  and  4  body  lengths 
behind  the  base  of  the  projectile.  The  grid 


4 


points  in  the  normal  direction  were  exponetially 
stretched  away  from  the  surface  with  the  minimum 
spacing  at  the  wall  of  .000020.  This  spacing 
locates  at  least  two  points  within  the  laminar 

sublayer. 

Th.  grid  shown  in  Figure  6  was  generated  in 
two  segments.  First,  the  grid  in  the  outer  region 
is  obtained  using  an  elliptic  solver9  for  the 
ogive  portion  and  straight-line  rays  for  the  re¬ 
maining  portion  which  runs  all  the  way  to  down¬ 
stream  boundary.  Second,  the  grid  in  the  base 
region  is  obtained  simply  by  extending  the 

straight  lines  perpendicular  to  line  AC  down  to 
the  center  line  of  symmetry  (line  BO).  An  expon¬ 
ential  stretching  with  the  minimum  spacing  of 
.000020  at  line  AC  is  used.  It  should  be  noted 

that  the  same  minimum  spacing  .000020  is  specified 
on  both  sides  of  the  cut  (line  AC)  thus  maintain¬ 
ing  a  smooth  variation  of  grid  across  the  cut. 
This  spacing  could,  of  course,  be  increased  down¬ 
stream  of  the  base.  The  number  of  grid  points 
above  and  below  line  AC  is  the  same  (50  points) 
which  means  that  an  adequate  number  of  points  are 
located  in  the  base  region.  As  can  be  seen  in 

Figure  6,  the  grid  points  are  clustered  near  the 
nose-cylinder  junction  and  at  the  projectile  base 
where  appreciable  changes  in  flow  variables  are 
expected. 

The  free  stream  Reynolds  number  for  the 
series  of  computations  was  fixed  at  4.5  *  10° 
based  on  the  total  model  length.  The  computations 
are  started  from  free  stream  conditions  and  march¬ 
ed  in  time  to  obtain  the  steady  state  solution. 
The  initial  calculation  was  made  for  M  =  0.9. 
Previous  converged  solutions  were  then  used  as 
starting  conditions  for  additional  Mach  number 
runs  to  achieve  faster  convergence. 

Figures  7  and  8  show  the  distribution  of  the 
surface  pressure  coefficient,  Cp  as  a  function  of 

axial  position,  x/D.  Figure  8  shows  the  overall 
view  whereas  Figure  7  shows  the  distribution  in 
the  near  wake  region  of  the  base.  The  distribu¬ 
tion  over  the  projectile  surface  itself  is  shown 
in  both  these  figures.  The  value  of  Cp  beyond  x/0 

=  6  is  the  value  of  pressure  coefficient  along  the 
cut  AC.  Both  these  figures  indicate  the  shock 
waves  near  the  nose-cylinder  junction  and  near  the 
blunt  base.  Although  not  shown  in  these  figures, 
the  pressure  along  the  base  remains  fairly 
constant  (within  *.005  variation  in  Cp  values). 

The  series  of  Figures  9,  10  and  11  show  the 
velocity  vector  field  in  the  base  region  for  M  = 
0.9  and  a  »  0.  Each  vector  shows  the  magnitude 
and  the  direction  of  the  velocity  at  that  point. 
Figure  9  shows  the  velocity  field  in  the  entire 
base  region.  One  can  see  the  expected  velocity 
defect  in  the  far  wake  region.  Figures  10  and  11 
show  the  velocity  field  in  the  vicinity  of  the 
base  (near  wake  region).  The  difference  between 
these  two  plots  being  that  the  former  one  is 
stretched  (not  drawn  to  same  scale)  while  the 
latter  is  drawn  with  the  same  scale  in  x  and  y 
directions.  Both  Figures  clearly  show  the  recir¬ 
culation  region  of  flow  in  the  base  reqion  and 
indirate  a  strong  shear  layer  as  well. 


to  the  same  scale  in  x  and  y  while  Figure  12  is 
no*.  nowever,  both  of  these  figures  are  drawn  to 
show  the  recirculation  region  and  the  position  of 
the  dividing  stream  line  as  clearly  as  possible. 
They  also  show  the  reattachment  point  which  for 
this  case  is  about  2  calibers  down  from  the  base. 

A  more  critical  check  of  the  computat i ona 1 
results  is  presented  in  Figure  14  where  the  base 
drag  is  plotted  as  a  function  of  Mach  number. 
Computational  results  are  indicated  by  circles, 
experimental  results7  by  triangles  and  the  squares 
indicate  the  results  obtained  using  a  semi- 
empirical  technique  developed  by  McCoy9.  Base 
drag,  as  expected,  increases  as  the  Mach  number 
increases  from  0.9  to  1.2.  The  sem> -empi r i ca 1 
technique  shows  generally  higher  base  drag  when 
compared  with  computational  and  experimental 

results.  The  computational  results  predict  the 

expected  drag  rise  that  occurs  for  0.9  <  M  <  1.2. 
The  computational  results,  however,  indicate  a 
greater  increase  in  drag  than  predicted  by  either 
the  semi-empirical  code  or  the  experimental  meas¬ 
urements.  The  discrepancy  between  the  numerical 
and  the  experimental  results  can  partly  be  attri¬ 
buted  to  the  fact  that  the  experimental  data  was 
obtained  with  a  sting  attached  to  the  base.  The 
sting  has  an  effect  of  weakening  the  recirculatory 
flow  in  the  base  region  and  leads  to  higher  base 
pressure  and  hence,  lower  base  drag. 

Figures  15,  16  and  17  show  the  variation  of 
pressure  drag,  skin  friction  drag  and  the  total 
drag  with  Mach  number  respectively.  The  rise  in 
the  pressure  drag  with  Mach  number  is  predicted 
correctly.  Skin  friction  drag  decreases  as  Mach 
number  i 'creases.  The  total  drag,  as  expected, 
increases  as  Mach  number  increases  from  0.9  to 
1.2.  The  computational  results  are  compared  with 
the  results  obtained  by  semi-empirical  technique 
and  the  agreement  is  cons.dered  good. 


V .  Summary 

A  procedure  has  been  described  in  which  the 
Azimuthal-Invariant  (generalized  ax i symmetric) 
thin-layer  Navier-Stokes  code  is  modified  in  such 
a  way  as  to  compute  the  base  flow  field  of  projec¬ 
tiles  at  transonic  speeds. 

The  computed  results  show  the  qualitative 
features  of  the  flow  field  in  the  base  region, 
namely  the  recirculation  region,  dividing  stream 
l>e,  reattachment  point,  etc.  Quantitative  com¬ 
parisons  of  the  base  drag  have  been  made  with 
other  available  data  for  various  Mach  numbers  in 
the  transonic  speed  range.  These  results  indicate 
that  the  present  numerical  technique  can  be  used 
successfully  to  predict  the  base  drag  of  projec¬ 
tiles  at  transonic  speeds. 

The  computed  results  for  this  paper  represent 
the  first  application  of  thin-layer  Navier-Stokes 
computational  technique  to  predict  projectile  base 
flow  at  transonic  velocity  using  the  flow  field 
segmentation  described  above.  The  results  indi¬ 
cate  that  this  technique  shows  good  promise  of 
providing  a  useful  new  computat ional  capability 
for  exterior  ballistics  of  shells. 


Future  computat i onal  efforts  will  investigate 
the  implementation  of  viscous  boundary  condition 


The  next  two  Figures  1?  and  13  are  stream 
function  contour  plots  in  the  near  wake  region, 
again  for  M  --  0.9  and  a  0.  Figure  13  is  drawn 


on  the  projectile  base,  improved  grid  resolution, 
and  alternate  turbulence  models. 
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Figure  1.  Spark  Shadowgraph  of  Secant-Ogive- 

Cylinder-Boattail  Projectile,  M  =  0.95, 

a  -  0.0 
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Figure  2.  Ax i symmetric  Body  and  Coordinate  System 
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Figure  3.  Schematic  Illustration  of  Flow  Field 
Segmentation 
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Figure  13.  Stream  Function  Contours,  M  =  0.9, 
a  *  0 


Figure  16.  Variation  of  Viscous  Drag  Coefficient 
with  Mach  Number,  a  *  0 
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Fiqure  IB.  Variation  of  Base  Drag  Coefficient 
with  Mach  Number,  a  =  0 


Figure  17.  Variation  of  Total  Drag  Coefficient 
wi th  Mach  Number ,0*0 
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Abstract 

A  computational  capability  has  been  developed 
for  predicting  the  flow  field  about  the  entire 
projectile,  including  the  recirculatory  base  flow, 
at  transonic  speeds.  Additionally,  the  computer 
code  allows  mass  injection  at  the  projectile  base 
and  hence  is  used  to  show  the  effects  of  base 
bleed  on  base  drag.  Computations  have  been  made 
for  a  secant-ogive-cylinder  projectile  for  a 
series  of  Mach  numbers  in  the  transonic  flow 
regime.  Computed  results  show  the  qualitative  and 
quantitative  nature  of  base  flow  with  and  without 
base  bleed.  The  reduction  in  base  drag  with  base 
bleed  is  clearly  predicted  for  various  mass  injec¬ 
tion  rates  and  for  Mach  numbers  .9  <  M  <  1.2.  The 
encouraging  results  obtained  indicate  that  this 
computational  technique  may  provide  useful  design 
guidance  for  shells  with  base  bleed. 

Nomenclature 

a  speed  of  sound 

A  cross  sectional  area  at  the  base 

Aj  injection  area  for  base  bleed 

CDb  base  drag  coefficient,  2  0b/p„u*A 

cp  specific  heat  at  constant  pressure 

C^  pressure  coefficient,  2(p  -  pj/p.u* 

D  body  diameter  (57.15mm) 

Dp  base  drag 

e  total  energy  per  unit  volume/ p^a^ 

E,  F,  q  flux  vector  of  transformed  Navier-Stokes 
.  equations 

H  n-invariant  source  vector 

I  mass  injection  parameter,  mj/pjJ.A 

J  Jacobian  of  transformation 

m,  mass  flow  rate  for  air  injection  at  the 
base,  PjUjA. 

M  Mach  number 

p  pressure/p„a* 

Pr  Prandtl  number,  w  C  /« 

R  body  radius  p 

Re  Reynolds  number, 

S  viscous  flux  vector 

t  physical  time 

u.v.w  Cartesian  velocity  components/a, 

U.V.W  Contravariant  velocity  components/aw 
x,y,z  physical  Cartesian  coordinates 


a  angle  of  attack 

y  ratio  of  specific  heats 

x  coefficient  of  thermal  conductivity 

u  coefficient  of  viscosity 

C,n,c  transformed  coordinates  in  axial, 

circumferential  and  radial  directions 
p  density/ p., 

t  transformed  time 

4  circumferential  angle 

Superscript 

*  critical  value 

Subscript 

b  base 

j  jet  conditions 

J  longitudinal  direction 

L  normal  direction 

o  total  conditions 

st  stagnation  conditions 

-  free  stream  conditions 


I.  Introduction 

A  major  area  of  concern  in  shell  design  is  in 
the  total  aerodynamic  drag.  The  designer,  ever 
desirous  of  increasing  the  range  and/or  terminal 
velocity  of  projectiles,  is  eager  to  decrease  the 
aerodynamic  drag. 

The  total  drag  for  projectiles  can  be  divided 
into  three  components:  (1)  pressure  (wave)  drag; 

(2)  viscous  (skin  friction)  drag;  and  (3)  base 
drag.  For  a  typical  shell  at  M  =  .90  the  relative 
magnitudes  of  the  aerodynamic  drag  components  are: 
(1)  pressure  drag,  20%;  (2)  viscous  drag,  30%,  and 

(3)  base  drag,  50%.  The  pressure  and  viscous  com¬ 
ponents  generally  cannot  be  reduced  significantly 
without  adversely  affecting  the  stability  of 
shell.  Recent  attempts  to  reduce  the  total  drag 
are  therefore  directed  at  reducing  the  base  drag. 

A  number  of  studies  have  been  made  to  examine 
the  drag  reduction  due  to  the  addition  of  a  boat- 
tail.  Although  this  is  very  effective  in  reducing 
the  drag,  it  has  a  negative  impact  on  the  aerody¬ 
namic  stability  of  shell  especially  at  transonic 
velocities.  An  excellent  review  of  base  drag  and 
the  effect  of  boattailing  is  presented  in 
Reference  1. 
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Another  effective  means  of  reducing  the  base 
drag  Is  that  of  'base  bleed1  or  'base  Injection'. 
In  this  method  a  small  amount  of  mass  Is  Injected 
into  the  base  region  which  increases  the  base 
pressure  and  thus  reduces  the  base  drag.  Recent 
range  and  precision  tests2  of  a  155nm  projectile 
with  and  without  base  bleed  have  been  conducted 
and  an  85%  reduction  in  base  drag  was  obtained. 
Presently  the  XM864  is  an  active  projectile  design 
which  is  attempting  to  use  the  base  bleed  concept 
for  increased  range.  This  concept  of  mass  injec¬ 
tion  at  the  projectile  base  has  been  widely 
studied  for  supersonic  flows  and  much  of  the  work 
has  been  reported  in  Reference  3.  One  limited 
study  at  supersonic  speeds  was  made  at  BRL  and  the 
results  were  reported  by  Dickinson1*.  A  limited 
study  made  in  the  transonic  flow  regime  has  been 
reported  in  Reference  5.  The  supersonic  regime 
has  typically  been  the  area  where  increased  range 
due  to  drag  reduction  has  been  studied.  Thus, 
only  limited  attention  has  been  focussed  on  the 
'base  bleed'  problem  in  transonic  flow. 

Most  of  the  work  using  the  ‘base  bleed' 
concept  has  been  either  experimental  or  semi- 
empirical  in  nature.  Sophisticated  numerical 
techniques  have  not  yet  been  utilized  to  predict 
the  effects  of  base  bleed  on  the  base  drag  reduc¬ 
tion.  Limited  computational  work  has  been  report¬ 
ed  recently  by  Sullins,  et  al6.  Their  work  dealt 
with  the  numerical  computation  of  the  base  region 
flow  of  a  supersonic  combustion  ramjet  engine 
using  two-dimensional  Navier-Stokes  equations. 
They  computed  the  flow  field  in  the  vicinity  of 
the  base  with  parallel  gas  injection  and  estab¬ 
lished  the  effect  of  base  injection  on  such  flows. 
Because  of  the  recent  advances  in  computer  techno¬ 
logy,  numerical  computational  capabilities  have 
been  developed  to  predict  the  aerodynamic  behavior 
of  artillery  shells.  Recent  papers7’8  have 
reported  the  development  and  application  of  the 
Azimuthal- Invariant  Thin-Layer  Navier-Stokes  comp¬ 
utational  technique  to  predict  the  flow  about 
slender  bodies  of  revolution  at  transonic  speeds. 
This  technique  has  been  modified  for  base  flow 
analysis  and  the  resulting  new  numerical  capabi¬ 
lity*  Is  used  here  to  predict  the  base  pressure  of 
shell  at  transonic  speeds  including  the  effect  of 
base  bleed.  Computed  results  show  the  quantita¬ 
tive  and  qualitative  details  of  the  base  flow 
structure.  The  technique  used  computes  the  full 
flow  field  over  the  projectile  at  transonic 
speeds;  therefore,  all  three  components  of  the 
total  drag  (pressure,  viscous,  and  base  drag)  are 
computed.  This  computational  technique  Is  then 
applied  to  predict  the  effects  of  base  bleed  on 
the  base  drag  reduction  at  transonic  speeds. 

A  brief  description  of  the  physical  problem 
and  the  governing  equations  Is  given  In  Sections 
II  and  III.  The  computational  technique  and  the 
method  of  solution  are  discussed  In  Section  IV. 
In  Section  V  results  are  shown  for  transonic  base 
pressure  computations  for  a  6-callber  secant- 
ogive-cylinder  shape  for  .9  <  M  <  1.2  with  and 
without  base  bleed.  Velocity  vector  plots,  stream 
function  contour  plots  and  density  contour  plots 
are  presented  to  show  the  qualitative  features  of 
the  flow  field  in  the  base  region.  Quantitative 
comparisons  of  base  drag  and  the  total  drag  both 
with  and  without  base  Injection  have  been  made. 
The  encouraging  results  show  that  the  present 
computational  technique  can  be  used  to  study  the 
effects  of  base  bleed  on  base  drag  and  can  possi¬ 


bly  have  a  positive  impact  on  the  XM864  devleop- 
ment.  Although  results  here  are  presented  for 
transonic  speeds,  current  computational  efforts 
are  directed  at  supersonic  velocities. 


Physical  Problem 


The  physical  problem  deals  with  the  transonic 
flow  over  a  projectile  shape  including  the  base 
region.  Although  the  entire  projectile  flow  is 
computed,  the  emphasis  here  is  on  the  flow  field 
in  the  base  region  of  the  projectile.  A  small 
amount  of  air  is  injected  at  the  projectile  base 
in  the  direction  parallel  to  the  primary  flow. 
The  injection  at  the  base  can  be  concentrated  at 
the  center  of  the  base  or  spread  throughout  the 
entire  base.  In  the  present  work,  however,  the 
injection  takes  place  over  90%  of  the  base. 
Figure  1  shows  a  schematic  illustration  of  the 
base  region  flow  field  with  base  injection.  The 
dividing  streamline  separates  the  reclrculary  base 
flow  from  the  primary  external  flow.  The  flow 
field  is  dominated  by  separation  and  mixed  regions 
of  locally  supersonic  and  subsonic  flows. 


The  complete  set  of  time- dependent  generali¬ 
zed  axisymmetric  thin-layer  Navier-Stokes  equa¬ 
tions  is  solved  to  obtain  a  numerical  solution  to 
this  problem.  The  numerical  technique  used  is  an 
implicit  finite-difference  scheme.  Although  time- 
dependent  calculations  are  made,  the  transient 
flow  is  not  of  primary  interest  at  the  present 
time.  The  steady  flow  is  the  desired  result  which 
is  obtained  in  a  time  asymptotic  fashion. 


Governi ng  Equatl ons 


The  Azimuthal  Invariant  (or  Generalized  Axi¬ 
symmetric)  thin- layer  Navier-Stokes  equations  for 
general  spatial  coordinates  E,  n,  c  can  be  written 
as7 

»Tq  +  *  H  «  Re'^S  (1) 


where  £  *  C(x,y,z,t)  Is  the  longitudinal 
coordinate 

h  *  n(y,z,t)  is  the  circumferential 
coordinate 

C  *  t(x,y,z,t)  Is  the  near  normal 
coordl nate 

t  •  t  is  the  time 


The  vector  of  dependent  variables  q  and  the  flux 
vectors  E,  G,  H  are  given  as 
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=  J 
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pW 

PuW+  CxP 
PvW+4yP 
pwW+CzP 
(<*p)W-CtP 


0 

0 

H  =  J"1^  pV[R5(U-tt)  +  Rc(W-ct)] 
-PVR*n(V-nt)  -  p/(R«n) 

0 

The  thin  layer  viscous  terms  valid  for  high 
Reynolds  number  flow  are  contained  in  the  vector 

S,  where 

r  0  1 


projectiles  or  swirl  flow  to  be  accomplished. 
There  is  some  evidence  which  indicates  that  base 
pressure  can  change  due  to  the  spin  of  a  projec¬ 
tile.  Although  the  present  work  considers  base 
flow  with  no  spin,  base  flow  with  spin  is  of 
interest  and  can  be  studied  using  the  present 
technique. 

For  the  computation  of  turbulent  flows  a  tur¬ 
bulence  model  must  be  supplied.  In  the  present 
calculations  a  Cebeci-type  two  layer  algebraic 
eddy  viscosity  model  as  modified  by  Baldwin  and 
Lomax10  is  used.  In  their  two  layer  model  the 
inner  region  follows  the  Prandtl-Van  Driest  formu¬ 
lation.  Their  outer  formulation  can  be  used  in 
wakes  as  well  as  in  attached  and  separated  bound¬ 
ary  layers.  In  both  the  inner  and  outer  formula¬ 
tions  the  distribution  of  vorticity  is  used  to 
determine  length  scales  thereby  avoiding  the 
necessity  of  finding  the  outer  edge  of  the  bound¬ 
ary  layer  (or  wake).  The  magnitude  of  the  local 
vorticity  for  the  axisymmetric  formulation  is 
given  by 


“(vy{;\  +  (x/3)(  txV'yV'A1'* 
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In  determining  the  outer  length  scale  a 
function10 


F(y)  *  y |“l  [1  -  exp(-y+/A+)]  (5) 


is  used  where  y+  and  A+  are  the  conventional 
boundary  layer  terms.  For  the  base  flow  (or  wake 
flow)  the  exponential  term  of  Equation  (5)  is  set 
equal  to  zero. 


The  velocities 

U  -  +  Cxu  +  tyv  ♦  ezw 

V  *  nt  *  V  +  V  +  nzw  <z) 

W  -  ct  +  Cxu  ♦  cyv  +  4zw 


represent  the  contravariant  velocity  components. 

The  Cartesian  velocity  components  (u,v,w)  are 
nondimensionalized  with  respect  to  a  (the  free 
stream  speed  of  sound).  The  density  (p)  Is  refer¬ 
enced  to  pm  and  total  energy  (e)  to  p„a2.  The 
local  pressure  is  determined  using  the  equation  of 
state, 

P  *  (Y-l)[e  -  0.5p(u** v*+w2)]  (3) 


where  r  1$  the  ratio  of  specific  heats. 

In  Equation  (1)  the  thin-layer  approximation 
Is  used  and  the  restrictions  for  axisymmetric  flow 
are  Imposed.  The  details  can  be  found  In 
References  8  and  9  and  are  not  discussed  here. 
Equation  (1)  contains  only  two  spatial  deriva¬ 
tives;  however  It  retains  all  three  momentum  equa¬ 
tions  thus  allowing  a  degree  of  generality  over 
the  standard  axisymmetric  equations.  In  parti¬ 
cular,  the  circumferential  velocity  Is  not  assumed 
to  be  zero  thus  allowing  computations  for  spinning 


IV.  Numerical  Method 
a.  Computational  Algorithm 

An  implicit  approximate  factorization 

finite-difference  scheme  in  delta  form  is  used  as 
described  by  Beam  and  Warming11.  An  Implicit 
method  was  chosen  because  it  permits  a  time  step 
much  greater  than  that  allowed  by  explicit 

schemes.  For  problems  in  which  the  transient 

solution  Is  of  no  interest,  this  offers  the  pos¬ 

sible  advantage  of  being  able  to  reach  the  steady 
state  solution  faster  than  existing  explicit 
schemes. 


The  Beam-Warming  implicit  algorithm  has 
been  used  in  various  applications7"13.  The  algo¬ 
rithm  can  be  first  or  second  order  accurate  In 
time  and  second  or  fourth  order  accurate  in  space. 
The  equations  are  factored  (spatially  split)  which 
reduces  the  solution  process  to  one-dimensional 
problems  at  a  given  time  level.  Central  differ¬ 
ence  operators  are  employed  and  the  algorithm  pro¬ 
duces  blodt  tridiagonal  systems  for  each  space 
coordinate.  The  main  computational  work  Is  con¬ 
tained  In  the  solution  of  these  block  trldiagonal 
systems  of  equations. 

b.  Finite  Difference  Equations 

The  resulting  finite  difference  equa¬ 
tions,  written  In  delta'form  are 


A. 
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( 1+ hfi^n-e  ,J‘ )  ( 1+ h«;Cn-e  jj‘ 

-  We'l4j.J'1HnJ)  *  (qrHl-qn)  »  -At{  S^A^G0  (6) 
-Re'l41Sn)-AtHn.eEJ-l[(S5Ai)2  ♦  (V?A;)2]Jqn 


Here  h  =  At  because  only  first  order  accuracy  in 
the  time  differencing  is  needed  for  the  steady 
state  flows  which  are  considered  here.  This 

choice  corresponds  to  the  Euler  implicit  time  dif¬ 
ferencing.  The  6’s  represent  central  difference 
operators,  A  and  V  are  forward  and  backward  dif¬ 
ference  operators  respectively.  The  Jacobian 

*  3E  9G 

matrices  A  =  — C  -  — x  along  with  the  coeffi- 

3q  3q 

cient  matrix  M  obtained  from  the  local  time 

linearization  of  $  are  described  in  detail  in 
Reference  6.  Fourth  order  explicit  (e^)  and 

implicit  (tj)  numerical  dissipation  terms  are 

incorporated  into  the  differencing  scheme  to  damp 
high  frequency  growth  and  thus  to  control  the  non¬ 
linear  instabilities.  A  typical  range  for  the 
smoothing  coefficients  is  eF  =  (1  to  5)  At  with 


c.  Flow  Field  Segmentation 

Figure  2  is  a  schematic  illustration  of 
the  flow  field  segmentation  that  is  used  to 
compute  the  entire  projectile  flow  field  including 
the  base  flow.  It  shows  the  transformation  of  the 
physical  domain  into  the  computational  domain  and 
the  details  of  the  flow  field  segmentation 
procedure  in  both  the  domains. 

The  cross  hatched  region  represents  the 
projectile.  The  line  BC  is  the  projectile  base 
and  the  region  ABCD  is  the  base  region  or  the 
wake.  The  line  AB  is  a  computational  cut  through 
the  physical  wake  region  which  acts  as  a  repeat- 
itive  boundary  in  the  computatioal  domain. 
Implicit  Integration  is  carried  out  in  both  E  and 
C  directions  (see  Figure  2).  Note  the  presence  of 
the  lines  BC  (the  base)  and  EF  (nose  axis)  in  the 
computational  domain.  They  both,  however,  act  as 
boundaries  in  the  computational  domain  and  special 
care  must  be  taken  in  inverting  the  block  tridiag- 
onal  matrix  in  the  E  direction.  The  details  are 
presented  in  the  next  section. 

d.  Implementation  of  Boundary  Conditions 
1.  Base  Flow  Without  Base  Injection 

The  no  slip  boundary  conditions  for 
viscous  flow  is  enforced  by  setting 


U  *  V  *  U  =  0  (7) 


on  the  projectile  surface  except  for  the  base.  At 
the  projectile  base  the  velocity  component  normal 
to  the  base  is  set  to  zero,  i.e.  U  »  0,  while 
other  flow  variables  are  set  equal  to  those  at 
grid  point  next  to  the  base.  In  other  words,  slip 


is  allowed  along  the  base  (inviscid  boundary 
condition).  Future  work  will  be  directed  at  the 
implementation  of  viscous  boundary  condition  at 
the  base  to  further  access  this  approximation. 

Along  the  computational  cut  (line  AB) 
the  flow  variables  above  and  below  the  cut  were 
simply  averaged  to  determine  the  boundary  condi¬ 
tions  on  the  cut.  On  the  centerline  of  the  wake 
region,  a  symmetry  condition  is  imposed. 


3u 

3Z 
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=  0 


(8) 


Free  stream  conditions  are  used  at 
the  outer  boundary.  Simple  extrapolation  for  all 
flow  variables  is  used  at  the  downstream  boundary 
(lines  AD  and  AG).  During  transient  calculations, 
use  of  a  specified  outflow  pressure  can  give  rise 
to  numerical  oscillations  in  the  base  region  flow 
field.  Eventually,  these  grow  and  swamp  the  solu¬ 
tion.  This  difficulty  is  avoided  by  simply  extra¬ 
polating  pressure  to  the  downstream  boundary  which 
is  the  procedure  always  used  with  supersonic  out¬ 
flow.  A  combination  of  extrapolation  and  symmetry 
is  used  at  on  the  nose  axis  (line  EF). 

As  a  result  of  the  flow  field  segmen¬ 
tation  procedure  described  in  Section  IV  c,  the 
block  tridiagonal  matrix  in  the  E  direction  has 
elements  at  J  -  JB,  JB+1  which  are  treated  as 
internal  boundaries  in  the  computational  domain  (J 
»  JB  represents  the  projectile  base  and  J  =  JB*1 
is  the  nose  axis).  The  block  tridiagonal  matrix 
in  the  E  direction  takes  the  following  form  (after 
setting  Ej  =  0  to  simplify  the  Illustration) 
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Here  A's  denote  the  quantity  A  and  I  is  a  SxS 

identity  matrix.  Note  the  appearance  of  two 
uncoupled  block  tridiagonals.  The  rows  at  JB  and 
JfHl  are  particularly  simple  because  boundary  con¬ 
ditions  are  updated  explicitly  at  the  end  of 
inversions.  These  changes  were  easily  implemented 
in  a  modular  fashion  into  an  existing  code  for 
projectile  base  flow  computations.  One  simply 
fills  the  block  tridiagonal  matrix  ignoring  the 
base  JB  and  the  nose  axis  JB+1.  Elements  in  these 
rows  are  then  overloaded  as  shown  above.  The  flow 
field  segmentation  does  not  affect  the  block  tri¬ 
diagonal  matrix  in  the  c  direction. 

2.  Base  Flow  with  Base  Injection 

The  boundary  conditions  used  for  base 
flow  with  mass  addition  are  presented  here.  The 
boundary  conditions  along  the  projectile  surface, 
at  the  cut  and  downstream  boundary  all  remain  the 
same  as  previously  described.  Along  the  base 
boundary  the  following  boundary  conditions  are 
Imposed. 


'  s  *j  5  0 

w  -  wj  s  wjb_i  (grid  point  next  to  the  base) 
p  *  pj  *  pst 

The  stagnation  density  is  obtained  from  the 
following  relation. 


The  grid  shown  in  Figure  4  was  generated 
in  two  segments.  First,  the  grid  in  the  outer 
region  is  obtained  using  an  elliptic  solver14  for 
the  ogive  portion  and  straight-line  rays  for  the 
remaining  portion  which  runs  all  the  way  to  down¬ 
stream  boundary.  Second,  the  grid  in  the  base 
region  Is  obtained  simply  by  extending  the 
straight  lines  perpendicular  to  line  AB  down  to 
the  center  line  of  symmetry  (line  CD).  An  expon¬ 
ential  stretching  with  the  minimum  spacing  of 
,000020  at  line  AB  is  used.  It  should  be  noted 
that  the  same  minimum  spacing  .000020  is  specified 
on  both  sides  of  the  cut  thus  maintaining  a  smooth 
variation  of  grid  across  the  cut.  This  spacing 
could,  of  course,  be  increased  downstream  of  the 
base.  The  number  of  grid  points  above  and  below 
line  AB  is  the  same  (B0  points)  which  means  that 
an  adequate  number  of  points  are  located  in  the 
base  region.  As  can  be  seen  in  Figure  4,  the  grid 
points  are  clustered  near  the  nose-cylinder  junc¬ 
tion  and  at  the  projectile  base  where  appreciable 
changes  in  flow  variables  are  expected. 

As  indicated  in  Figure  4,  the  fine  viscous 
grid  follows  the  cut  labeled  as  AB  in  Figure  2. 
In  so  far  as  the  viscous  shear  layer  begins  to 
neck-down  shortly  behind  the  base,  much  of  this 
fine  grid  resolution  is  wasted.  As  a  consequence 
logic  has  been  implemented  to  adjust  the  grid  cut 
AB  to  the  viscous  shear  layer.  Such  a  grid  is 
shown  in  Figure  5  in  which  the  height  of  the  cut 
is  determined  from  a  moment  of  shear  subject  to 
various  constraints  and  averaging.  Specifically, 

the  cut  height,  "zj  at  each  J-location  is  deter¬ 
mined  by  the  relation 
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The  amount  of  air  injected  into  the  base  region 
can  be  specified  by  the  mass  flow  rate,  m^.  Since 
Pj  and  Aj  are  known,  uj  can  be  calculated  for  any 
given  mass  flow  rate.  Rather  than  specifying  irij 
however,  it  is  customary  to  specify  a  mass  injec¬ 
tion  parameter,  I  where  I  *  mj/pjj^A. 

e.  Computational  Grid 

The  finite  difference  grid  used  for  the 
numerical  computations  was  obtained  from  a  grid 
generator  developed  in  Reference  14.  This  program 
allows  arbitrary  grid  point  clustering,  thus  enab¬ 
ling  grid  points  for  the  projectile  shapes  to  be 
clustered  In  the  vicinity  of  the  body  surface. 
The  grid  consists  of  108  points  In  the  longitudin¬ 
al  direction  and  50  points  In  the  radial  direc¬ 
tion.  The  full  grid  Is  shown  In  Figure  3  while 
Figure  4  shows  an  expanded  view  of  the  grid  in  the 
vicinity  of  the  projectile.  The  computational 
domain  extended  to  4  body  lengths  In  front,  4  body 
lengths  In  the  radial  direction  and  4  body  lengths 
behind  the  base  of  the  projectile.  The  grid 
points  in  the  normal  direction  where  exponentially 
stretched  away  from  the  surface  with  the  minimum 
spacing  at  the  wall  of  .000020.  This  spacing 
locates  at  least  two  points  within  the  laminar 
sublayer. 


where  the  summation  is  carried  out  only  for  those 
points  within  an  interval  .2D  <  ZJL  <  0/2.  Here  D 

is  the  base  diameter,  S?  is  a  central  difference 

operator  and  e  is  a  positive  parameter  which 
ensures  a  standard  grid  i  f  all  4z  uJL  are  zero  or 

if  e  is  very  large.  Additional  averaging  is  used 
in  the  x-direction  (longitudinal  direction). 
Preliminary  results  have  been  obtained  using  the 
grid  shown  in  Figure  5  and  further  computations 
are  underway. 


V.  Results 

The  model  geometry  used  in  the  present  study 
is  shown  in  Figure  6.  The  model  consists  of  a  3 
caliber  secant-ogive  nose  and  a  3  caliber 
cylinder. 

The  free  stream  Reynolds  number  for  the 
series  of  computations  was  fixed  at  4.5  *  106 
based  on  the  total  model  length.  The  computations 
are  started  from  free  stream  conditions  and  march¬ 
ed  in  time  to  obtain  the  steady  state  solution. 
The  initial  calculation  was  made  for  M  *  0.9. 
Previous  converged  solutions  were  then  used  as 
starting  conditions  for  additional  Mach  number 
runs  to  achieve  faster  convergence.  The  results 
are  now  presented  for  *both  cases,  (1)  base  flow 
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without  Injection  end  (11)  base  flow  with 
injection. 

Figures  7  and  8  show  the  distribution  of  the 
surface  pressure  coefficient.  Cp  as  a  function  of 

axial  position  respecti vely ,  without  and  with  mass 

injection  at  the  base.  The  value  of  Cp  beyond  X/0 

=  6  is  the  value  of  pressure  coefficient  along  the 
line  extending  from  the  cylinder  portion  straight 
to  the  downstream  boundary.  In  Figure  7  there  is 
no  mass  injection  at  the  base.  The  pressure  dis¬ 
tribution  reflects  the  shock  pattern  that  typical¬ 
ly  occurs  on  shell  at  transonic  velocities,  the 
rapid  expansion  which  occur  at  the  blunt  base  and 
the  recompression  that  occurs  downstream  of  the 
base.  Although  not  shown  in  these  figures,  the 
pressure  along  the  base  remains  fairly  constant 
(within  t.005  variation  in  Cp  values).  The  pres¬ 
sure  coefficient  distribution  for  a  case  with 
large  mass  addition  is  shown  in  Figure  8.  The 
previous  rapid  expansion  at  the  base  and  recom¬ 
pression  downstream  of  the  base  are  seen  to  be 
virtually  eliminated. 

Figure  9  shows  the  velocity  vector  field  in 
the  base  region  for  M  =  0.9,  a  =  0  and  1  »  0. 
Each  vector  shows  the  magnitude  and  the  direction 
of  the  velocity  at  that  point.  The  figure  shows 
the  velocity  field  when  there  is  no  base  bleed  and 
the  recirculatory  flow  in  the  base  region  is 
clearly  evident. 

The  velocity  vector  plots  in  Figures  10,  11 
and  12  show  the  effect  of  base  bleed  on  the  near 
wake  flow  field.  Figure  10  shows  the  effect  of 
base  bleed  for  the  case  when  the  mass  injection 
parameter  is  very  small  (I  =  .01).  The  change  in 
the  flow  field  is  not  very  dramatic.  In  Figure  11 
the  mass  injection  parameter  is  increased  to  .07, 
and  the  effect  of  base  bleed  can  be  clearly  seen. 
The  near  wake  flow  field  has  changed  considerably. 
Figure  12  shows  the  effect  of  base  bleed  for  a 
still  higher  mass  injection  parameter,  I  =■  .13. 
The  flow  field  in  the  base  region  has  now  been 
dramatically  altered.  The  recirculation  pattern 
has  been  totally  swept  downstream. 

The  next  four  Figures  13,  14,  15  and  16  are 
stream  function  contour  plots  in  the  wake  region, 
again  for  M  =  0.9  and  o  *  0.  All  these  figures 
are  deliberately  stretched  in  y  direction  (not 
drawn  to  the  same  scale  in  x  and  y)  to  show  the 
flow  pattern  in  the  base  region  as  clearly  as  pos¬ 
sible.  Figure  13  is  for  the  case  of  base  flow 
with  no  mass  injection  at  the  base.  It  clearly 
shows  the  recirculation  region  and  the  position  of 
the  dividing  streamline  which  separates  the  recir¬ 
culatory  base  flow  from  the  main  flow.  The  reat¬ 
tachment  point  is  about  2  calibers  down  from  the 
base.  Note  the  strong  shear  layer  in  the  base 
region. 

Figures  14,  15  and  16  show  the  flow  pattern 
In  the  base  region  with  mass  Injection  allowed  at 
the  base.  Figure  14  shows  the  effect  of  base 
bleed  when  the  mass  injection  parameter  is  very 
small  (I  =>  .01).  The  reattachment  point  remains 
at  about  the  same  place  as  with  no  injection  at 
the  base.  The  flow  pattern  has  changed  slightly 
as  can  be  seen  by  the  dividing  streamline,  however 
the  recirculation  region  has  not  changed  dramatic¬ 
ally.  In  Figure  15,  the  mass  injection  parameter. 


1,  has  been  increased  to  .07  and  now  the  effect  of 
mass  injection  can  be  clearly  seen.  The  reattach¬ 
ment  point  has  moved  further  down  stream.  The 
flow  pattern  in  the  near  wake  flow  field  has 
changed  considerably  and  the  separation  bubble  is 
reduced  in  size.  When  the  mass  injection  para¬ 
meter  is  increased  further,  I  =  .13,  its  effect  on 
the  flow  field  in  the  base  region  is  apparent. 
Figure  16  shows  that  dramatic  change  in  the  flow 
field.  It  indicates  no  presence  of  any  recircula¬ 
tion  region  and  shows  how  the  shear  layer  has  been 
displaced  markedly. 

A  more  critical  look  at  the  computational 
results  is  presented  in  Figures  17  through  20. 
These  figures  show  the  quantitative  details  of 
projectile  flow  field.  Figure  17  shows  the  varia¬ 
tion  of  base  drag  with  mass  injection  rates  for  M 
=  0.9  and  a  =  0.  The  reduction  in  base  drag  with 
base  injection  can  be  seen  clearly.  The  percent 
reduction  in  base  drag  increases  with  the  an 
increase  in  the  injection  rate. 

Since  the  entire  projectile  flow  field, 
including  the  base  flow,  has  been  computed,  all 
three  drag  components  have  been  computed  and  thus 
the  total  drag  determined.  Figure  18  shows  the 
variation  of  the  total  drag  with  varying  mass 
injection  rates.  Again,  the  reduction  in  the 
total  drag  is  apparent.  As  the  injection  rate  is 
increased,  the  percent  reduction  in  total  drag 
increases. 

Figures  19  and  20  show  respectively,  the 
variation  of  base  drag  and  the  total  drag  with 
Mach  number  both  with  and  without  base  injection. 

In  both  these  figures  the  computational  results 
without  injection  at  the  base  are  shown  by  the 

solid  line  whereas  the  dotted  line  represent  the 

computational  results  obtained  with  injection  at 
the  base.  The  reduction  in  base  drag  and  also  the 
total  drag  with  base  injection  can  be  clearly 
seen.  Figure  19  indicates  that  the  percent  reduc¬ 
tion  in  base  drag  has  increased  with  an  increase 
in  Mach  number  from  .9  to  .98.  In  both  the 
figures  the  expected  drag  rise  in  the  transonic 

speed  regime  is  well  predicted  for  .9  <  M  <  1.2 

and  the  reduction  in  base  drag  and  the  total  drag, 
due  to  base  bleed  has  been  clearly  demonstrated. 


V.  Summary 

A  promising  computational  capability  has  been 
developed  which  computes  the  full  projectile  flow 
field  including  the  recirculatory  base  flow  at 
transonic  speeds  both  with  and  without  base 
injection. 

Numerical  computations  have  been  made  for 
Mach  numbers  .9  <  M  <  1.2  to  predict  the  base  drag 
and  the  total  drag  with  and  without  base  bleed. 
Computed  results  show  the  qualitative  features  of 
the  flow  field  in  the  near  wake  for  both  cases. 
The  effect  of  base  injection  on  the  qualitative 
nature  of  base  flow  has  been  clearly  shown.  Quan¬ 
titative  comparisons  of  base  drag  and  the  total 
drag  both  with  and  without  base  injection  have 
been  made.  For  M  1  0.9  and  a  =  0  the  computation¬ 
al  results  show  the  reduction  in  base  drag  and  the 
total  drag  for  several  mass  injection  parameters. 
Results  are  also  presented  for  .9  <  M  <  1.2  for  a 

given  mass  injection  rate  and  the  reduction  in 
* 


6 


base  drag  and  the  total  drag  has  been  demonstrated 
for  this  range  of  transonic  speeds. 


Current  computational  efforts  are  directed  at 
the  numerical  computation  of  base  flow  at  super¬ 
sonic  speeds.  The  encouraging  results  obtained 
thus  far  at  transonic  speeds  indicate  that  the 
computational  technique  presented  here  shows  the 
promise  of  providing  the  capability  to  predict  the 
base  drag  and  hence  the  total  drag  both  with  and 
without  base  injection. 
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Figure  1.  Schmatic  Illustration  of  Base  Region 
Flow  Field  with  Base  Bleed 
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Figure  13.  Stream  Function  Contours,  M  =  0.9 

a  =  0,  1=0 


Figure  16.  Stream  Function  Contours,  M  =  0.9, 
a  =  0,  I  =  .13 
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Figure  19.  Variation  of  Base  Drag  Coefficient 
with  Mach  number,  o  =  0  (with  and 
without  Base  Bleed) 
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Figure  20.  Variation  of  Total  Orag  Coefficient 
with  Mach  Number,  a  *  0  (with  and 
without  Base  Bleed) 
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III.  GRIDS 


(;<><><l  solution  accuracy  depends  on  having  properly  spaced  grids  that  are 
smoothly  varying  and  not  overly  skewed.  During  the  research  program  considerable 
ellorl  was  therefore  devoted  to  the  task  of  grid  generation.  A  projectile  grid  genera¬ 
tion  program  that  can  use  either  elliptic  or  hyperbolic  grid  generation  procedures 
was  devised.  The  attached  paper,  presented  at  the  1082  Army  Numerical  Analysis 
and  C  omputers  Conference,  describes  our  basic  grid  generation  solver. 
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ABSTRACT,  The  determination  of  accurate  projectile  aerodynamics  is  a 
major  area  of  concern  for  shell  designers  involved  with  new  shapes  and 
Ballisticans  concerned  with  developing  artillery  aiming  data.  To  achieve  the 
desired  goals  a  research  effort  has  been  on  going  within  the  Aerodynamics 
Research  Branch/BRL  to  establish  a  predictive  capability  for  determing  pro¬ 
jectile  aerodynamics.  Modern  finite  difference  codes  have  been  applied  to  the 
projectile  problem  and  encouraging  results  have  been  obtained  in  transonic1’2 
and  supersonic3  flow.  The  generation  of  good  computational  grids  has  been  a 
prerequisite  for  achieving  these  flow  field  solutions. 

This  paper  describes  a  versatile  grid  generation  program  which  has  been 
developed  for  standard,  hollow  and  non- ax i symmetric  projectile  shapes.  The 
grid  generator  makes  use  of  both  elliptic  and  hyperbolic  type  partial  differ¬ 
ential  equation  solvers.  The  code  allows  arbitrary  grid  point  clustering 

along  the  body  suface  in  areas  of  anticipated  flow  field  gradients.  The  outer 
boundary  can  also  be  arbitral ly  defined  with  its  own  clustering  distribution. 
The  grid  is  then  generated  between  these  two  boundarys  with  either  straight 
rays  or  by  use  of  an  elliptic  solver.  For  those  cases  when  the  outer  boundary 
is  not  restricted,  the  grid  can  be  generated  using  a  hyperbolic  solver  which 
adds  the  additional  benefit  of  an  orthogonal  mesh. 

The  mathematical  development  of  the  clustering  functions  and  partial 
differential  equation  solvers  are  described  and  a  series  of  grids  are  pre¬ 
sented  which  show  the  versatility  of  the  grid  generation  program.  Grids  for 
ogive-cylinder-boattai  1  configurations,  hollow  ring  airfoil  projectiles  and 
non-axi symmetric  projectiles  are  discussed. 

1.  INTRODUCTION.  The  numerical  solution  of  the  Navier-Stokes4 ,b 
equations  has  been  successfully  applied  to  a  wide  variety  of  problems.  The 
versatility  of  these  methods  is  inpart  attributed  to  the  solution  of  the 
transformed  set  of  differential  equations.  Using  transformed  equations  the 
physical  space  can  be  mapped  onto  a  regularly  spaced  rectangular  region  for 
two  dimensional  flow.  This  mapping  allows  for  a  wide  variety  of  projectile 
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configurations  to  be  solved  using  the  same  basic  numerical  technique.  An 
example  of  some  characteristic  projectile  shapes  are  shown  in  Figure  1.  A 
standard  projectile  shape  which  consists  of  an  ogive  cylinder  boattail  is 
shown  in  la;  a  more  non-conventional  shape  but  one  of  considerable  interest, 
the  triangular  boattail  configuration  in  lb;  and  a  tubular  projectile 
configuration  which  has  been  type  classified  and  is  currently  in  the  Army 
inventory,  in  lc.  To  calculate  the  flow  field  for  any  one  of  these  shapes  the 
first  requirement  is  to  develop  a  suitable  finite  difference  grid  for  use  with 
the  equation  solver.  The  grid  generator  described  in  this  paper  addresses 
this  problem. 

Grid  generation  routines  are  employed  to  generate  a  network  of  constant  Z 
and  n  lines  in  the  physical  x-y  plane  as  indicated  in  Figure  2a.  Correspond¬ 
ing  uniform  values  of  Z  and  n  in  the  computational  space  define  a  one  to  one 
mapping  between  points  j,k  in  the  physical  plane  to  points  j,k  in  the  computa¬ 
tional  plane  as  shown  in  Figure  2b.  The  mapping  functions  are  described,  at 
least  numerically,  once  ■,  and  are  known  in  the  physical  plane  as  a 

function  of  x  j  >(c  and  yj^.  The  metric  quantities  Z%,  Zy,  nx,  and  needed  in 

the  transformed  flow  equations  can  then  be  determined  numerically  (see,  for 
example,  References  4-6). 

The  grid  generation  program  presented  here  describes  earlier  work  done  by 
the  authors7  as  well  as  extensions  which  include  a  hyperbolic  solver  and  the 
addition  of  more  general  projectile  shapes.  The  grid  generator  is  modular  and 
begins  with  a  determination  of  the  body  shape.  The  inner  body  clustering 
routine  is  then  called  to  distribute  points  in  the  vicinity  of  previously 
determined  flow  field  gradients.  The  next  option  allows  for  the  insertion  of 
stings  for  wake  modeling,  a  rear  cut  or  forward  cut.  If  the  outer  boundary  is 
free  or  unconstrained  as  is  the  case  for  conventional  projectiles,  the  hyper¬ 
bolic  solver,  which  generates  a  smoothly  varying  orthogonal  grid,  is  called. 
For  those  cases  where  the  outer  boundary  is  constrained,  as  is  the  case  for 
tubular  projectile  shapes,  the  outer  boundary  clustering  routine  is  called. 
Once  the  outer  boundary  is  specified  the  elliptic  solver  is  called.  The  grids 
generated  up  to  this  point  would  be  planar  and  sufficient  for  axi-symmetric 
calculations.  However  for  three  dimensional  flow  fields  a  periodic  or  non¬ 
periodic  grid  is  generated  by  spinning  the  planar  grid  about  the  symmetry 
axis.  A  flow  chart  of  the  overall  grid  program  is  shown  in  Figure  3. 

The  following  sections  of  the  paper  will  present  some  of  the  details  used 
for  the  inner  boundary  clustering  the  outer  boundary  description  and  interior 
grid  generation. 

2.  INNER  BOUNDARY  DESCRIPTION.  The  body  shape  can  be  input  to  the 
program  by  cards,  file  specification  or  as  a  set  of  x,y  ordinates.  The  data 
is  assumed  to  be  non-dimensional  with  respect  to  the  diameter  or  cord  depend¬ 
ing  on  the  projectile  configuration.  Additionally,  the  code  can  generate  a 
parabolic  arc  or  standard  class  of  projectiles  such  as  sharp  or  blunt,  tangent 
or  secant  ogive-nose,  cylindrical  body,  boattail,  or  spherical  cap.  Once  the 
body  shape  is  determined  the  values  of  x  along  the  body  axis  are  distributed 
by  contiguously  combining  segments  of  the  clustering  function 


Xj  .  X0  *  a»j  4  b*j  4 
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where  ^  =  ( j-jQ)/( jf-jQ)  and  j  is  an  index  value  such  that  points  jQ  to  jf 
lie  in  the  interval  xQ  to  Xf  and  xj  =  xQ  while  xj^.  =  Xf.  Equation  (1)  is 

used  to  cluster  Xj  as  a  function  of  j.  The  user  determines  the  shape  of  the 

clustering  function  by  specifying  the  initial  and  final  increments  of  x,  that 
is 
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Since  xQ  and  Xf  are  also  specified,  a,  b,  and  c  are  determined 

c  =  {Vxf  +  AxQ  -  2h(xf  -  xQ) }/( h  -  3h2  +  2h3) 
b  =  (Axq  -  h(xf  -  xQ)  -  c(h3  -  h) }/( h2  -  h) 
a  =  xf  -  xQ  -  b  -  c 

where  h  =  (jf  -  jQ)_1. 

The  amount  of  clustering  at  each  point  is  determined  by  the  specified 
values  of  axq  and  Vxf.  Moreover,  because  Axq  and  Vx^  are  specified,  the  user 

can  smoothly  patch  functions  together  to  form  a  general  clustering  function. 
One  drawback  to  the  clustering  function,  Eq.  (1),  is  that  the  function  is  not 
guaranteed  to  be  monotone  in  the  interval.  This  can  happen,  for  example,  if 
axq  is  too  small  and  Vxf  too  large. 

At  this  point  a  sting  or  forward  cut  can  be  added  to  the  previously 
described  body  as  shown  in  Figures  4a  and  4b.  Again  the  clustering  function 
of  Equation  (1)  is  used  to  distribute  points  along  these  new  boundaries. 

3.  GRID  GENERATION  USING  A  HYPERBOLIC  SOLVER.  For  most  projectile 
applications  the  outer  boundary  is  unconstrained  and  simply  needs  to  be  placed 
far  enough  away  from  the  projectile  body  so  as  not  to  adversely  affect  the 
flow  field  solution.  This  situation  represents  an  ideal  case  for  a  hyperbolic 
grid  generation  scheme. 
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Once  the  body  points  have  been  redistributed  and  the  sting  or  cut  has 
been  determined,  a  grid  can  be  generated  using  a  hyperbolic  solver  similar  to 
that  described  In  Reference  8.  Before  the  actual  solver  can  be  Implimented 
however,  the  distance  to  the  outer  boundary  must  be  specified  and  either  con¬ 
stant  spacing  in  n  or  some  type  of  stretching  function  is  required.  The  n 
stretching  used  here  is  determined  by  the  following  relationship 


')k-l.k  ■  l.«W-  1  <3> 

Here  asq  is  the  minimum  specified  grid  spacing  desired  at  the  wall  or  inner 

boundary.  The  parameter  e  is  determined  by  a  Newton-Raphson  iteration  process 
so  that  the  sum  of  the  above  increments  matches  the  known  arc  length  between 
n  =  0  and  n  =  nmax  for  points  which  have  the  same  value  of  £. 

The  governing  equations  for  the  hyperbolic  solver  are  obtained  by 
requiring:  (1)  the  coordinate  lines  C  and  n  to  be  orthogonal;  and  (2)  the 
specification  of  a  cell  volume  or  area  for  the  two  dimensional  case.  The 
condition  of  orthogonality  requires 


AC  •  An  =  0 


(4) 


The  second  equation  is  obtained  by  specifying  a  grid  cell  volume  (or  area  in 
two  dimensions).  Since  the  grid  cell  volume  is  finite  the  transformation 
Jacobian  will  be  greater  than  one,  i.e.. 


dxdy  =  |x^yn  -  x^J  dCdn  (5) 

The  set  of  grid  generation  equations  are  therefore  given  in  the  physical  plane 
by 


C  n  +  C  n  =0 
x  'x  y  y 

C  n  -  5  n  =  J 
x  y  y  x 

or  in  the  transformed  plane  by  (6) 

Vn  *  y5yn  -  0 

Vn  '  V«  *  lyJ  5  V 

Using  local  linearization  for  this  set  of  non-linear  differential  equations, 
the  resulting  system  is  shown  to  be  hyperbolic8  and  can  therefore  be  marched 
in  the  n  direction. 
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The  linearized  set  of  differential  equations  to  be  solved  numerically  is 
written  in  vector  form  as 


Ar  +  Br*  =  ? 
4  n 


(7) 


where 


LV*V“J  '  LyJ 

where  x°,  y°,  etc.,  refers  to  known  conditions. 

The  set  of  Equations  (7)  are  solved  with  an  implicit  finite  difference 
scheme  which  is  first  order  accurate  in  the  n  direction(k)  and  where  central 
differencing  is  used  in  the  4  direction(  j) .  The  resulting  set  of  finite 
difference  equations  becomes 


A(rj+l,k+l  '  rj-l  ,k+l^  +  B(rj,k+1,  ~  rj,k)  =  *j,k+l  (8) 

2A4  An 

Rearranging  Eq.  (8)  and  setting  An  =  A4  =  1  results  in 

1  rj+l  ,k+l  +  8  rj,k+l  "Z  rj-l  ,k+l  =  ?j,k+l  +  8  rj,k  =  3j,k+l  (9) 


where 


j  .k+1 


(x°x°  +  y°y°) 


(_yoxo  + 


x0y0)j>k  +  v  +  V° 


Equation  (9)  is  now  in  a  form  which  can  be  easily  solved  by  inverting  a  block 
tridiagonal  matrix  with  2  *  2  blocks.  The  terms  x°  and  y°  are  central  differ¬ 
enced  as  4  4 
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,0  _  Xj+l,k  "  Xj-l,k 


*j.k 


(10) 


.  yj+l,k  ~  yj-l,k 


j  ,k 


The  terms  x°  and  y°  are  obtained  from  Equation  (6)  evaluated  at  the  old 
n  n 


station(o).  That  is 


xoxo  +  yoyo  =  o 


X0V0  _  x 0  y 0  =  wO 


(ID 


Solving  for  x°  and  y°  with  x°  and  y°  given  in  (10)  yields 


-  v° 
«“  ■ — i—7 


n  (x°.y°) 


x?  V° 

y°  = - 1 - r 

n  (x°+y°) 


(12) 


1 


The  cell  volume  remains  to  be  specified.  This  specification  is  important 
since  it  has  the  effect  of  controlling  the  grid  evolution  as  the  solution  is 
being  marched  out  from  the  body.  The  method  chosen  here  is  straight  forward 
and  uses  the  stretching  function  given  by  Equation  (3).  Specifying  the 
minimum  spacing  at  the  wall  AsQ  and  the  total  number  of  points,  jmax,  in 

the  n  direction  an  array  of  arc  lengths  Ask  is  determined.  Since  the  Ax  is 

known  along  the  j  line,  the  volumes  are  calculated  by 


V  =  (Ask)  (x.j+1>k  -  *j>k)  (13) 

This  specification  of  cell  volumes  yields  smoothly  varying  grids  in  the 
n  direction.  Grid  volume  control  is  obtained  by  varying  the  arc  length 
distribution  Ask  and/or  surface  point  distribution.  An  additional  volume 

specification  approach  can  be  found  in  Reference  8.  A  grid  generated  using 
this  technique  is  shown  in  Figure  5a  and  5b  for  a  standard  projectile 
configuration  with  sting. 

4.  OUTER  BOUNDARY  DEFINITION.  For  those  cases  where  the  outer  boundary 
is  constrained  or  specified  a  grid  point  distribution  along  the  outer  boundary 
is  required.  An  example  is  shown  in  Figure  6  .  A  part  of  the  grid  generation 
problem  then  is  the  formation  of  an  arbitrary  outer  boundary.  Here  this 
boundary  is  built  up  by  connecting  contiguous  cubic  segments,  which  in  the 
degenerate  case  can  be  straight  lines.  Figures  7a  and  7b  illustrate  two 
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typical  outer  boundary  curves.  In  Figure  7 a  three  cubic  segments  make  up  the 
boundary,  n  =  n__„.  Each  segment  is  formed  by  specifying  the  values  of  x,y, 

and  angle  6,  at  the  endpoints,  where  0  is  the  angle  between  the  curve  and  the 
x  axis.  In  the  example,  Figure  7a,  0  =  90°  ,  0.  =  0  =  0° ,  or  180°  and  o  = 


The  data  (x,y,0)  at 
curves 


which  are  equivalent  to 


each  endpoint  determines 

ctjt  + 

Sjt  +  62t2 
cubi  c 


the  shape  of  the  parametric 

o  <  t  <  1  (14) 


y  =  yQ  +  Vx-Xq)  +  Vx'x0)2  +  y3(x*x0)3  U5) 

The  parametric  cubic  is  used  because  the  condition-^  =  <*>  can  be  specified 
(segment  be  of  Figure  7b  has  this  constraint  at  both  endpoints). 

The  solution  for  the  parameters  a,,  ,  8.,  and  B0  can  be  found  in  Refer¬ 
ence  7.  1212 

The  outer  boundary  curve  is  thus  made  up  of  contiguous  cubic  segments 
starting  from  the  5=0  boundary.  Points  are  distributed  along  this  curve 
either  as  a  uniform  distribution  of  arc  length,  or  as  a  specified  arc  length 
distribution  using  the  previously  defined  clustering  scheme,  Eq.  (1).  Since 
the  true  arc  length  is  not  specified  a  priori,  precise  alignment  of  points 
along  the  outer  boundary  can  be  determined  only  after  the  cubic  segments  are 
specified  and  the  arc  length  is  computed. 

5.  STRAIGHT  RAY  AND  ELLIPTIC  GRID  GENERATION.  Once  the  boundary  curves 
have  been  specified  and  points  are  distributed  on  the  n  =  0  and  n  bound¬ 
aries,  two  types  of  grid  generation  procedures  can  be  used.  n,a* 

In  the  first  case,  lines  of  constant  5  (i.e.,  the  rays  emerging  from  the 
body)  are  formed  by  simply  connecting  straight  lines  from  points  along  n  =  0 
to  points  along  n  =  ’Vlax*  ^he  sPac’n9  '•n  n  a1°n9  each  such  line  is  either 

uniform  or  is  determined  by  the  stretching  relationship  given  by  Equation 
(3).  Figures  8a  and  8b  illustrate  a  straight  ray  grid  with  clustering  in  n 
for  a  tubular  projectile. 

In  the  second  case,  the  grid  is  generated  with  elliptic  partial  differen¬ 
tial  equations  following  References  9,  10,  and  11.  The  grid  generating  equa¬ 
tions  are  solved  on  the  specified  computational  space  for  unknowns  x,  l  and 
yj  .k  • 
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(16) 


ax,r  -  2fix.  +  yx  =  -J2  (Px.  +  "0x  ) 
££  £n  nn  '  £  ir 

o  _ 

aycr  -  20yr  +  yy  =  -J  ("Pyr  +  Qy  ) 
J££  J£n  Jnn  v  J 

where 


2  2 
=  x  +  y  , 
n  n 


0  = 


x  rx  +  y  ry  , 
£  n  J  £J  n 


Y  = 


2  2 
X£  +  y£  ’ 


J  = 


x  fy 

v  n 


x  yr 

n-7  £ 


and 


•p  =  p  +  P  e“a^  n_nmax^ 

o  m 

g  =,  gQ  e-»("-n0)  ,  g^  .-bln-n^) 


Here  PQ,  Q0,  Pm,  Qm,  a  and  b  are  prescribed  clustering  parameters.  Along  the 
n  =  0  and  n  =  n„  w  boundaries,  x;  t,  and  y^  >.  have  been  previously  prescribed. 

nia X  J  r  J  j* 

Along  the  £  =  0  and  £  =  E  .  which  are  either  vertical  or  horizontal  lines  in 

iTldX 

the  physical  space,  the  following  boundary  conditions  are  enforced:  either 

x  is  given  and  -  0 


on  a  vertical  boundary,  or 


(17) 


x^  =  0  and  y  is  given 


on  a  horizontal  boundary. 

The  difference  equations  to  Eq.  (16)  {see  Reference  7)  are  solved  with  a 
successive  line  over  relaxation  (SLOR)  procedure.  As  an  initial  guess  for  the 
relaxation  procedure  the  straight  line  ray  procedure  previously  described  is 
used.  For  the  most  part,  if  coefficients  P  and  Q  are  large,  the  SLOR  pro¬ 
cedure  is  very  difficult  to  converge.  Consequently,  the  algebraic  clustering 
function,  Eq.  (3)  is  recommended. 

In  the  algebrajc  clustering  approach  the  elliptic  solver  is  used  to  gen¬ 
erate  a  grid  with  P  =  Q  =0.  The  x,y  points  along  a  £  =  constant  line  are 
then  redistributed  along  this  line  as  a  function  of  arc  length.  The  clus¬ 
tering  function  Eq.  (3)  is  used  for  this  purpose.  This  procedure  works  quite 
well  and  provides  excellent  control  of  the  grid  spacing  near  the  body  surface. 
Further  details  are  given  in  Reference  7. 
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The  elliptic  solver  need  not  be  used  over  the  entire  range  in  Because 
of  the  boundary  condition,  Eq.  (17),  the  elliptic  equations  can  be  joined  to  a 
straight  ray  along  any  vertical  or  horizontal  boundary  line  in  5.  Figure  9 
shows  details  of  such  a  procedure  used  for  a  secant-ogive-cylinder  boattai  1 
projectile  which  also  includes  a  sting.  Here  the  C-region  over  the  secant- 

ogive  nose  is  generated  using  the  elliptic  equations  while  the  remainder  is 
meshed  with  straight  rays.  After  the  basic  grid  is  formed,  the  entire  grid  is 
clustered  in  n  using  Eq.  (3). 

6.  3D  GRIDS.  The  final  option  available  in  the  code  is  the  ability  to 

generate  three  dimensional  grids.  At  present  the  grids  are  formed  in  a  two 
dimensional  plane  and  then  rotated  about  a  symmetry  axis.  The  rotation  is 

either  periodic  or  non-periodic  depending  on  the  grid  desired.  For  cases 
where  the  flow  field  has  planar  symetry,  such  as  a  projectile  at  anqle  of 

attack,  without  spin,  a  non-periodic  grid  is  generated. 

The  generation  of  grids  for  projectile  shapes,  with  non-axi symmetri c 
sections  (Figure  lb)  is  accomplished  with  a  series  of  planar  grids.  Planes 
are  generated  normal  to  the  projectile  axis  at  incremental  values  of  Ax.  For 
each  of  these  planes  a  grid  is  generated  using  an  0  type  grid  (Figure  10). 
These  grids  are  then  combined  to  form  a  three  dimensional  mesh  making  sure 

that  continuity  in  the  x  direction  is  maintained. 

7.  SUMMARY.  A  versatile  grid  generation  program  has  been  described 
which  utilizes  general  elliptic  and  hyperbolic  equation  solvers  for  internal 
grid  generation.  The  flexibility  of  longitudinal  grid  point  distribution  is 
obtained  with  the  general  clustering  functions  allowing  points  to  be  placed  in 
the  vicinity  of  flow  field  gradients.  Grid  clustering  is  also  obtained  near 
the  body  surface  for  viscous  flow  field  calculations. 

A  series  of  grids  have  been  presented  which  show  the  versatility  of  the 
code.  Grids  for  secant-ogi ve-cylinder  boattails  have  been  shown  using  an 
elliptic  solver,  hyperbolic  solver  and  a  hybrid  elliptic/straight  ray  solver. 
The  generation  of  a  grid  for  a  non-conventi onal  hollow  projectile  shape  has 
been  demonstrated. 
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a.  Conventional  Secant-Ogive  Cylinder  Boattail  (SOCBT)  Projectile 


b.  Secant-Ogive  Triangular  Boattail  Projectile 
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c.  Tubular  Projectile 


Figure  1.  Projectile  Configurations 
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C  =  C  (x,y) 
17s17  («,y ) 


(x.y  SPECIFIED) 


vl0  (x,y  SPECIFIED) 


figure  2.  Mapping  from  Physical  Space  to  Computational  Space 
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Figure  3.  Flow  Chart  of  Grid  Generator  Programs 
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Figure  5b.  Expanded  View 


Figure  6.  Tubular  Projectile  Grid  or  C-Grld 
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Figure  7.  Outer  Boundary  Structure  and  Terminology 
for  Two  Classes  of  Grid 
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Figure  9a.  Hybrid  Elliptic  and  Straight  Ray  Grid  for  SOCBT  With  Sting 
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Figure  10a.  Projectile  Cross  Section  with  Periodic  B.C.  (0-Grid) 
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Figure  10b.  Projectile  Cross  Section  with  Symmetry  Plane  (0-Grid) 
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IV.  NUMERICAL  ALGORITHMS 


A  majority  of  tin*  research  ellort  was  devoted  to  numerical  algorithm  work  so  as 
to  improve  computational  efficiency.  Much  of  this  work  is  described  in  two  invited 
survey  papers  which  are  attached. 


-  1  - 


Implicit  Finite  Difference  Simulation 
of  Inviscid  and  Viscous  Compressible  Flow 


presented  at  the  Symposium  on 
Transonic,  Shock  and  Multidimensional  Flows 
May  13-15,  1981,  Madison,  Wisconsin 


IMPLICIT  FINITE  DIFFERENCE  SIMULATION 


Of 

INVISCID  AND  VISCOUS  COMPRESSIBLE  FLOW 

JOSEPH  L.  STEGER 


I.  INTRODUCTION. 

It  is  not  always  convenient  to  use  the  simplified  equa¬ 
tions  that  extract  the  essential  physics  from  the  more  com¬ 
plete  set  of  inviscid  and  viscous  fluid  conservation- law-equa¬ 
tions.  Such  a  situation  may  occur  if  the  usually  inviscid 
outer  flow  is  highly  rotational  and/or  if  the  viscous  layer  is 
fully  separated. 

Numerical  procedures  for  solving  the  system  of  conserva¬ 
tion-law-equations  of  fluid  flow  are  not  as  efficient  as,  say, 
the  numerical  procedures  developed  for  the  scalar  nonlinear 
potential  equation  used  in  inviscid  transonic  flow  analysis. 

Of  course,  the  solution  of  a  system  of  equations  requires 
more  work  than  the  solution  of  a  scalar  equation.  Equally  or 
more  significant,  however,  is  that,  in  dealing  with  a  system 
of  equations  one  often  encounters  characteristic  speeds  (i.e. 
eigenvalues)  of  disparate  magnitude  (stiffness)  and  of  both 
positive  and  negative  sign.  These  last  conditions  can  make 
use  of  implicit  differencing  schemes  desireable  and  put  rath¬ 
er  severe  constraints  on  the  choice  of  spatial  differencing 
operators. 

The  purpose  of  this  paper  is  to  review  the  use  of  implic¬ 
it  finite  difference  schemes  to  solve  the  Euler  and  Navier- 
Stokes  equations  in  primitive  variables.  In  part  one  of  this 
paper  an  approximate  factorization  (AF)  implicit  finite  dif¬ 
ference  scheme  for  solving  the  Euler  and  Navier-Stokes  equa- 


tions  is  discussed.  The  equations  are  cast  in  generalized  co¬ 
ordinates  and  partial  differential  equation  grid  generation 
techniques  are  used.  In  this  approach  the  flux  vectors  of  the 
equations  are  differenced  as  whole  quantities  and  time-accu¬ 
rate  or  time-like  iterative  schemes  are  used  to  solve  the 
equations  for  general  boundary  surfaces.  In  part  two  of  this 
paper  ways  of  splitting  and  reducing  the  governing  equations 
are  reviewed  with  an  aim  towards  developing  more  accurate,  ef¬ 
ficient,  and  robust  numerical  algorithms.  Again,  implicit 
schemes  are  emphasized.  Here,  though,  the  methods  are  less 
developed. 

II.  IMPLICIT  FINITE  DIFFERENCE  FLOW  FIELD  SIMULATION. 

Over  the  last  several  years,  a  set  [1-6]  of  versatile, 
somewhat  robust  computer  codes  has  been  developed  for  simulat¬ 
ing  steady  or  unsteady  inviscid  or  viscous  compressible  flow. 
The  computer  programs  make  use  of  general  coordinate  transfor¬ 
mations,  numerical  grid  generation  techniques,  viscous  model¬ 
ing,  and  implicit  finite  difference  algorithms  to  achieve  a 
high  degree  of  adaptiveness  to  flow  conditions.  In  this  sec¬ 
tion  a  review  of  this  overall  methodology  is  put  forth.  For 
brevity,  the  discussion  here  is  restricted  to  two-dimensional 
compressible  flow,  although  the  basic  procedures  have  also 
been  applied  to  three-dimensional  flow  [2],  incompressible 
flow  [7],  and  supersonic  flow  solved  by  parabolic-like  march¬ 
ing  [  8  ] . 

a)  Transformed  Thin-Layer  Equations 

As  governing  equations  [1,9]  we  take  the  two-dimensional 
thin-layer  Navier-Stokes  equations  subject  to  general  coordi¬ 
nate  transformation,  but  kept  in  conservation-law-form  [10,11] 

3  Q  +  3-F  +  3  G  »  Re_13  S  (2.1) 

t  t  h  n 


where 


£  =  £ (x,y,t) 


n  =  n (x,y,t) 


t  =  t 


and  the  flux  terms  are  defined  as 
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where 


mx  ■  (2/3)u(2nxun  -  nyvn) 

m2  *  w(Vn  +  nxV 

m3  *  <2/3)u(2nyvn  -  nxun) 

m4  -  <P~l(Y  -  D'^x^Ce2) 
®5  -  <P;X(Y  -  l)“1ny3n(c2) 


Here  p  is  density,  p  is  pressure,  u  and  v  are  Cartesian  veloc¬ 
ity  components,  and  c  is  the  sound  speed.  The  total  energy 
per  unit  volume,  e,  is  defined  by 

e  -  (y  -  l)_1p  +  0 . 5p (u2  +  v2)  (2.2) 


and  the  (unsealed)  contravariant  velocities  are  defined  as 

U  =  5t  +  ?xu  +  ?yv  (2.3a) 

v  =  nt  +  nxu  +  nyv  (2.3b) 

The  metrics  5^.,  etc. ,  are  determined  once  a  mapping 

is  defined.  Usually,  a  numerical  mapping  is  employed.  The 
metrics  are  related  to  xT ,  x^ ,  etc.,  by  the  relations 


Kx  ”  JV 

5y  -  -JV 

5t  *  "Vx  - 

nx  -  -Jy5 

'  hy  =  JX^ , 

nt  =  "xtnx  '  yrn; 

-i 

J  *  xry 

~  xnyC 

Here  €  varies  around  the  body  surface,  and  n  varies  away  from 
the  body  surface,  as  indicated  in  Fig.  1-3.  The  symbol  ~  de¬ 
notes  that  the  scalar  or  vector  quantity  is  divided  through  by 
the  Jacobian,  J. 

For  practical  viscous  flow  calculations,  a  turbulence 
model  is  needed.  The  algebraic  two-layer  eddy-viscosity  model 
as  proposed  by  Baldwin  and  Lomax  [9]  is  used.  The  thin-layer 
approximation  requires  that  Re  >>  1  and  that  the  body  coincide 
with  an  n  »  const  line. 

The  inviscid  part  of  the  governing  equations  is  kept  in 
conservation  law  (i.e.  divergence)  form  so  as  to  capture  as 
accurately  as  possible  the  Rankine  Hugoniot  shock  jump  rela¬ 
tions.  Conservation-law-form  is  also  useful  in  implicit  for¬ 
mulations  in  that  it  can  lead  to  cleaner  local  linearization 
formula.  Conversely  it  can  cause  numerical  inaccuracy  unless 
metric  transformation  terms  are  properly  dealt  with. 

b)  Comments  About  the  Transformed  Equations 

The  transformed  equations  offer  several  significant  ad¬ 
vantages  over  the  less  complicated  Cartesian  form  of  the  equa¬ 
tions.  Chief  among  these  is  that  fact  that  the  physical 


boundary  surfaces  can  coincide  with  transformed  coordinate 
lines.  This  feature  can  be  used  to  simplify  the  application 
of  boundary  conditions.  Body  conforming  coordinates  are  also 
necessary  to  simplify  the  governing  equations  so  as  to  permit 
the  use  of  the  thin-layer  viscous  model.  Another  significant 
aspect  of  coordinate  transf ormations  is  that  they  can  be  used 
to  cluster  grid  points  to  flow  field  action  regions. 

Generally  we  prefer  to  use  a  single  transformation  of  co¬ 
ordinates  to  map  the  physical  plane  onto  a  uniform  rectangular 
computational  plane.  Ideally  the  body  boundary  surfaces 
should  lie  on  the  boundaries  of  the  computational  domain.  In 
this  way  one  arrives  at  a  well-ordered  system  of  finite  dif¬ 
ference  equations  which  is  sparse,  which  can  be  efficiently 
solved  for,  and  which  is  usually  amenable  to  vectorized  com¬ 
puter  processing.  These  advantages  are  so  significant  that 
they  are  relinquished  with  great  sorrow.  The  use  of  trans¬ 
forms  and  finite  difference  methods  in  no  way  requires  such 
a  simple  topology,  however.  Although  seldom  used,  it  is  pos¬ 
sible  to  match  or  overlap  more  than  one  grid  or  coordinate 
system  to  treat  complex  geometries  or  to  remove  grid  related 
stiffness  of  the  equations,  sketches  1  and  2  illustrate.  In 
using  such  multiple  grid  systems,  one  will  have  to  deal  with  a 
more  complex  program  and  take  care  not  to  introduce  numerical 
instability  at  grid  interfaces. 


c)  Metric  Accuracy 

Although  the  conservation- law- form  of  the  equations  is 
useful  for  capturing  shocks  and  helps  simplify  the  local  lin¬ 
earization  process,  its  use  can  lead  to  inaccuracy  in  certain 
difference  formulations.  In  putting  (2.1)  into  conservation- 
law-form,  use  is  made  of  the  exact  relation 


Ft^yn  -  3ny?]  +  Gt"85xn  +  3nx5I  *  0 


(2.5) 


In  differencing  the  flow  equations,  (2.1),  operators  5^  and 

<5^  are  introduced  to  approximate  3^  and  3^.  If  for  example, 

the  metrics  y  ■  £  /J  etc.  are  exactly  evaluated,  then  (2.5) 
n  X 

cannot  be  equal  to  zero  but  is  equal  to  a  small  error,  the 


truncation  error.  This  error  shows  up  in  (2.1)  as  a  source  or 
sink,  a  "metric  tare".  For  a  highly  stretched  grid  the  metric 
tare  in  differencing  (2.1)  can  be  so  appreciable  as  to  invali¬ 
date  the  approximation  and  even  numerical  instability  can  re¬ 
sult. 

To  avoid  this  problem  steps  must  be  taken  as  described  in 
[1,2,5].  As  first  noted  in  [1],  if  £  /J  =  y  etc.  are  cen- 

X  T] 

trally  differenced  with  the  same  operators  5,  and  6  used  to 
difference  the  fluid  flux  terms,  then  the  difference  equations 
also  exactly  satisfy  (2.5).  Central  differencing  the  metrics 
removes  the  source  term  difficulty,  although  flow  field  solu¬ 
tion  accuracy  can  still  be  poor  if  a  very  poor  grid  is  used. 
This  process  can  be  extended  to  three-dimensions  [2,5]. 

As  an  alternative  to  the  above,  good  success  has  been  ob¬ 
tained  by  making  what  we  have  termed  "free  stream  subtraction" 
[2].  In  this  approach  (2.1)  is  put  into  a  perturbation-like 
form 

3tQ  +  3?(F  -  FJ  +  3n(G  -  GJ  -  Re'^S  (2.6) 

That  is,  the  "metric  tare"  is  approximately  subtracted  off. 
This  works  on  a  well  generated  grid  because  the  metric  tare  is 
usually  only  severe  in  the  far  field. 

If  one  is  willing  to  accept  a  "weak  conservation- law  - 
form",  the  equations  could  be  modified  as  (here  for  inviscid 
flow  only) 

AAA  A 

3XQ  +  3?F  +  3^G  -  FO^-3^)  +  G  (-S^x^+S^x^)  =  H  (2.7) 

/v 

The  right  hand  side  source  term  H  should  not  effect  the  shock 
strength  or  location,  but  it  will  contribute  diagonal  terms 
which  by  themselves  can  be  weakly  unstable. 

d)  Grid  Generation 

To  take  advantage  of  the  transformed  governing  equations 
it  is  necessary  to  generate  a  smoothly-varying  body-conforming 
grid.  While  this  is  a  difficult  task  in  general,  a  variety 
of  algebraic  [c.f.  12-14]  and  partial  differential  equation 


[c.f.  14-22]  schemes  have  been  developed  which  at  least  handle 
the  two-dimensional  grid  generation  about  a  single  body  in  a 
fairly  automatic  way.  By  automatic  is  meant  that  if  the  user 
carefully  specifies  grid  points  and  clustering  information 
along  the  mesh  boundaries,  the  grid  generation  schemes  will 
usually  generate  a  smooth  nonsingular  interior  grid.  Of 
course  the  grid  may  not  be  as  optimum  as  one  would  like,  but 
the  results  are  likely  to  be  satisfactory.  Moreover,  progress 
is  being  made  in  treating  more  complex  two  and  three-dimen¬ 
sional  configurations. 

The  grid  generation  methods  based  on  solving  partial  dif¬ 
ferential  equations  have  always  seemed  especially  appealing. 
This  is  because  the  numerical  expertise  one  develops  for  sol¬ 
ving  the  flow  equations  is  directly  applicable  to  the  grid 
generation  task.  Figures  la  and  lb  show  an  example  of  a  grid 
using  elliptic  partial  differential  generating  equations  as 
taken  from  [23],  while  Figs.  2a  and  2b  show  a  grid  obtained 
from  hyperbolic  partial  differential  generating  equations  as 
taken  from  [22].  In  both  cases  the  grid  lines  are  orthogonal 
to  the  body  and  the  grid  spacing  at  the  body  is  uniformly  con¬ 
trolled.  The  grid  generated  with  hyperbolic  partial  differen¬ 
tial  equations  is  ideal  for  many  external  flow  configurations, 
it  is  essentially  orthogonal  throughout. 

e)  Difference  Equations  and  Numerical  Algorithm 

An  implicit,  noniterative,  time-accurate  finite  differ¬ 
ence  algorithm  has  been  used  to  solve  the  transformed  govern¬ 
ing  equations.  Although  viscous  flows  are  ideally  treated 
with  an  implicit  scheme,  the  same  numerical  algorithm  is  used 
for  inviscid  flow  calculations  as  well.  By  doing  so  one  can 
base  the  time  step  size  on  accuracy  considerations  and  not  be 
overly  concerned  about  highly  clustered  or  distorted  meshes . 
Computer  programs  that  use  explicit  or  semi-implicit  (e.g.  ex¬ 
plicit  in  the  streamwise  direction,  [c.f.  24])  schemes  can  be 
more  efficient  for  a  given  problem,  but  are  generally  not  as 
versatile. 


The  Beam-Warming  delta  form  approximate  factorization  al¬ 
gorithm  (25,26],  with  various  adaptations,  has  been  used  to 
solve  the  thin-layer  equations.  It  is  remarked  that  similar 
numerical  algorithms  have  been  developed  independently,  and  in 
aerodynamics  applications  the  contributions  of  Briley  and 
MacDonald  (27,28]  are  notable.  For  either  trapezoidal  or 
Euler  implicit  temporal  differencing  the  delta  form  differen¬ 
cing  scheme  for  the  thin-layer  equation  is  given  by: 

(I  +  h«SrAn  -  J_1e.h7-AcJ)  (I  +  h<5_Bn  -  J-1e .  hV  A  J 
t  i  t  t  n  inn 

-Re_1hJ  J_1Mn)  (Qn+1  -  Qn)  -  -At(5,Fn  +  <5  Gn  -  Re-1?  Sn) 
n  t  n  n 


-eehJ~1[ (V^A^) 2  +  (VnAn)2]JQn 


(2.8) 


Here  h  *  X+£*  a  *  0  or  1  for  first  or  second  order  accuracy, 
and  £.,  e  *  0(1)  with  e-  >  2e  are  added  numerical  dissipa- 
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tion  terms.  The  operators  6^,  6^  are  three  point  central  dif¬ 
ferences  e.g. 


Vi  -  Fi-1 


while  7,  A  are  conventional  backward  or  forward  operators. 


V  -  Qj  -  Qj-i* 


Finally,  T  is  the  midpoint  operator 


Sk+*i  "  Sk-» 


where  S  itself  contains  midpoint  differences  so  the  viscous 

A  A  A 

term  uses  three  points.  The  coefficient  matrices  A,  B,  and  M 


result  from  local  linearization  of  the  flux  terms,  with  A  and 

A  A  A  A  A  A 

B  the  Jacobian  matrices  [3F/3Q]  and  [3G/3Q]  while  M  contains 
derivative  operators  from  linearizing  S. 

The  difference  equations  given  by  (2.8)  are  readily 
solved  with  an  alternating  direction-like  sequence  in  which 
one  inverts  all  block  tridiagonals  in  5  followed  by  all  block 
tridiagonals  in  n.  Boundary  conditions  can  be  imposed  implic¬ 
itly,  approximate  implicitly,  or  explicitly  if  the  implicit 
stability  range  is  not  upset.  These  and  additional  details, 
including  use  of  higher  order  space  differencing,  are  de¬ 
scribed  in  [1-7,25,26], 

f )  Results 

The  numerical  scheme  described  previously  has  been  used 
in  a  variety  of  steady  and  unsteady  state  flow  problems  as  re¬ 
ported  in  [1-7,23],  The  simulation  of  aileron  buzz  represents 
a  typical  application  of  the  code,  so  a  few  such  calculations 
are  presented  below.  These  results  are  reproduced  from  [4], 

In  the  buzz  calculations  a  rigid  aileron  is  allowed  a  one 
degree  of  freedom  motion  about  its  hinge  line  as  described  by 
the  equation 

Ia6a  =  H(t)  (2'9) 

Here  I  is  the  aileron  mass  moment  of  inertia  about  the  hinge 

cl 

line,  H  is  the  aerodynamic  forcing  term,  and  6  is  the  aileron 
angle  of  deflection,  see  Fig.  3.  Equation  (2.9)  is  solved 
along  with  the  thin- layer  equations  on  a  C-type  grid  in  which 
the  grid  deforms  to  follow  the  aileron  motion. 

The  results  of  a  numerical  simulation  for  a  NACA  65-213 
a  *  0.5  airfoil  section  are  indicated  in  Figures  4  to  8.  Ex¬ 
perimental  data,  albeit  with  some  three-dimensional  effect,  is 
available  on  this  airfoil  from  the  tests  of  Erickspn  and 
Stephenson  [29]  in  which  they  mounted  the  wing  of  the  P-80 
from  the  sidewall  of  the  Ames  sixteen  foot  tunnel. 


U  .  1'  J 


According  to  the  experiment,  at  MOT  =  0.82  and  a  =  -1  deg, 
the  aileron  could  be  restrained  at  an  angle  near  zero,  and 
when  freed,  would  buzz.  In  the  numerical  simulation  the  ai¬ 
leron  was  initially  deflected  to  4  deg,  it  would,  on  being 
released,  buzz  as  indicated  in  Fig.  4.  The  computed  frequency 
of  22.2  Hz  is  in  good  agreement  with  the  experimental  value  of 
21.2  Hz.  However,  in  the  numerical  calculation,  the  aileron 
deflects  ±11.1  deg  about  the  angle  -1.1  deg,  while  the  corre¬ 
sponding  experimental  values  are  ±9.2  deg  about  -3  deg. 

In  the  numerical  calculations  at  a  slightly  higher  Mach 
number,  Mm  =  0.83,  the  aileron  does  go  into  a  buzz  cycle  when 
freed  from  a  zero  deflection  position.  An  essentially  steady- 
state  solution  was  used  as  initial  data.  The  build-up  of  ai¬ 
leron  deflection  angle  as  a  function  of  time  is  indicated  in 
Fig.  5.  After  four  cycles,  a  quasi-steady-state  is  reached 
and  the  aileron  oscillates  at  22.7  Hz. 

At  a  lower  Mach  number,  M  =0.79,  the  aileron  motion 

00 

damps  to  a  neutral  value  even  though  the  flap  was  initially 
deflected  4  deg.  Data  from  this  calculation  are  displayed  in 
Fig.  6. 

Several  frames  from  computer-generated  film  strips  show¬ 
ing  Mach  contours  are  shown  in  Figs.  7a-7c  at  selected  times 
for  a  =  -1  deg  and  M^  =  0.82.  Contour  levels  of  M  =  0.2,  0.4, 
0.6,  0.98,  1.0,  and  1.02  were  used  in  order  to  illustrate  both 
the  separated  flow  regions  and  sonic  lines. 

Finally,  for  the  aileron  held  fixed  at  a  higher  Mach  num¬ 
ber,  Mm  =  0.85,  we  find  that  the  viscous  flow  does  not  reach 
a  steady  state  but  buffets  at  a  frequency  of  about  26.6  Hz. 
This  is  indicated  by  the  unsteady  hinge  moment  coefficient,  CH 
shown  in  Fig.  8.  If  the  aileron  is  then  released,  it  no  long¬ 
er  oscillates  in  a  simple  sinusoidal  motion.  Viscous  effects 
appear  to  be  much  more  dominant  and  change  the  frequency  and 
amplitude  of  the  aileron  motion.  Similar  type  motion,  but  for 
a  different  airfoil,  has  been  observed  experimentally  [30]. 


III.  FLUX  SPLIT  SYSTEMS. 


In  a  series  of  recent  papers  [31,26,32]  splitting  of  the 
spatial  flux  terms  of  the  conservative  form  of  the  Euler  equa¬ 
tions  was  proposed.  The  flux  terms  are  usually  split  based  on 
the  positive  and  negative  eigenvalues  (characteristic  speeds) 
of  their  appropriate  Jacobian  matrices.  Nonconservative 
schemes  based  on  the  same  principle  have  previously  been  pro¬ 
posed  [33,34],  and  have  recently  been  developed  for  aerodyna¬ 
mic  applications  [35-37].  Related  older  schemes  have  been 
identified  and  others  are  under  extensive  development  [c.f. 

38] . 

Splitting  the  flux  vectors  (conservative  form)  or  coef¬ 
ficients  matrices  (nonconservative  form)  based  on  the  sign  of 
their  eigenvalues  allows  the  use  of  upwind  (either  backward  or 
forward)  spatial  differencing  schemes.  Without  use  of  such 
splitting  only  central  spatial  differencing  can  be  used  to 
approximate  the  Euler  equation  flux  derivatives,  except  of 
course,  for  those  spatial  directions  in  which  the  coordinate 
velocity  exceeds  the  sound  speed.  Upwind  differencing  schemes 
can  offer  some  advantages  over  central  differencing  insofar 
that  they  are  mere  dissipative,  in  some  instances  can  follow 
the  physics  better,  and  can  lead  to  new  implicit  approximate 
factorization  schemes.  It  is  this  latter  property  which  is 
the  subject  of  this  section. 

The  aim  of  the  conservative  form  plus  and  minus  flux 
vector  splitting  is  to  recast  the  inviscid  portion  of  the 
equations  into 

3 .  Q  +  3  F+  +  3  F-  +  3  G+  +  3  G~  =  Viscous  Part  (3.1) 

t  x  x  y  y 

where  Cartesian  coordinates  are  used  for  illustration,  i.e. 
(2.1)  with  £  =  1  =  n  and  E  =  0  =  n  •  The  Jacobian  matrices 

+  y  y  x 

3F  /3Q  and  3G  /3Q  are  constructed  to  have  positive  real  eigen¬ 
values  while  3F~/3Q  and  3G~/3Q  are  to  have  negative  real  ei- 
genvalues.  The  initial  development  of  F"  and  G  relied  on  the 
fact  that  the  Euler  equations  are  homogeneous  of  degree  one 
and  proceeded  as  (illustrated  for  F) 


F  =  AQ  A  =  3F/3Q 

=  SAS_1Q 

=  S(A+  +  A")S_1Q  with  2A±  =  (A±|A|)  (3.2) 

=  SA+S-1Q  +  SA-S_1Q 
=  F+  +  F- 


where  S  is  a  matrix  of  the  eigenvectors  of  A  while  A  is  a  di¬ 
agonal  matrix  of  its  eigenvalues.  This  approach  doesn't  work 
exactly  as  desired  because  the  crucial  eigenvalues  are  those 

of  the  Jacobian  matrices  3F+/SQ  and  3F  /3Q,  and  3F±/3Q  5* 

+  _  1  + 

SA”S  .  Nevertheless  the  eigenvalues  of  3F  /3Q  have  the  pro¬ 
per  signs  if  not  the  proper  magnitudes  [39]. 

A  general  formula  for  the  flux  vectors  F“,  G”  is  gxven 
by  [32]: 


_P 

2y 


2(y-l)A1  +  X3  +  X4 
2(y-l)riu  +  X3(u+ck1)  +  r4(u-ck1) 
2(y-1)A1v  +  A3(v+ck2)  +  ^4(v-ck2) 
(y-1)I1(u2+v2)  +  [ (u+ck^) 2  +  (v+ck2)2] 

+  -J  [  (u-ck^) 2  +  (v-ck2)2]  +  WJ;[ 


where  k^  and  k2  =  1  or 


0  for  F"  or  0  and  1  for  G“ 


(3.3) 


(3  -y)  (A.  +  A.)c2 

w  =  _ Z _ * 

II  2(Y  -  1) 


and  (in  the  present  application) 


For  example,  the  eigenvalues  i.  of  F  are  u,  u,  u+c,  u-c,  and 
F“  is  defined  from  5^  using 

2A^-  =  u  i  { u | 

21^”  =  u  +  c  ±  |u+c|  (3.5) 

2A^”=u-c±  | u-c | 

4  #  + 

Similar  relations  hold  for  G“  with  v,  v,  v+c,  v-c.  In  gener- 
alized  coordinates  F  as  given  by  (2.1)  has  eigenvalues  [40,1] 
U,  U,  U  ±  c/p  2  _  2  and  F“  is  derived  from  (3.3)  with 

2A1±  -  U  ±  |0|  U  =  Cxu  +  ?yv  +  Ct 

2A3±  =  U  +  c  ±  | U+c I 

2A4±  =  U  -  c  ±  | U-c | 

where  kj_  -  5X/(CX2  +  and  k2  =  5y/ (Cx2  +  Cy2) *5- 

A  previously  identified  [32]  difficulty  with  the  above 
±  + 

formulation  is  that  F  and  G-  have  discontinuous  derivates 
because  | A |  has  a  discontinuous  derivative.  As  discussed  in 
[32]  it  is  necessary  to  smooth  } A | ,  and  this  is  neatly  accom¬ 
plished  [41]  by  replacing  (3.4)  with 


where  e  is  small.  As  indicated  in  Fig.  9,  this  gives  a 
smooth  A ^7  variation  which  asymptotes  to  the  old  formulation. 
It  also  adds  numerical  dissipation  whenever  an  eigenvalue 
changes  sign.  In  numerical  tests  on  a  one-dimensional  tran¬ 
sonic  nozzle  the  use  of  the  new  formulation  (3.7)  in  place  of 
(3.4)  gives  a  smooth  sonic  line  result  that  was  not  previous¬ 
ly  obtained,  see  Fig.  10. 


Other  split  flux  vectors  have  been  proposed.  For  exam¬ 
ple,  Bram  van  Leer  [41]  has  suggested  the  form  (here  given  in 
one-dimension) 

u  <  c 

^pc  (M+l )  2/4  (3.8) 

pc2(M+l)2[ (y-1)M+2]/(4y)  and  F-  *  F  -  F+ 

Y2[(f,+)2/f,+]/[2(Y2-l) ] , 

'  / 

F+  =  F  and  F~  =  0 

This  splitting,  devised  from  different  arguments  than  (3.2), 
is  naturally  smooth  at  points  where  the  eigenvalues  change 
sign.  In  the  above  test  problem  use  of  the  van  Leer  flux 
vector  gives  the  pleasing  result  shown  in  Fig.  11. 

Because  upwind  differencing  schemes  generate  lower  or 
upper  triangular  matrices,  flux  split  implicit  algorithms  can 
be  devised  for  the  inviscid  equations  which  are  efficiently 
inverted.  For  example,  a  second  order  fully  implicit  dif¬ 
ferencing  of  (3.1)  is  obtained  using  three  point  upwind  dif¬ 
ferencing  in  space  and  in  time 

(5tQ  +  6xF+  +  <SxF"  +  5yG+  +  4yG”,jk1  =  °  (3‘9) 

where 

6^Qn+1  -  (3Qn+1  -  4Qn  +  Qn_1)/(2At) 

6xFj  *  (_3Qj  +  4Qj+l  "  Qj+2>/<2Ax>  etc* 

With  use  of  local  linearization  to  avoid  iterative  solution 
of  the  nonlinear  terms,  and  with  use  of  approximate  factori¬ 
zation  to  simplify  the  inversion  work,  a  delta  form  implicit 
differencing  of  (3.9)  can  be  obtained  as 


(I  +  h6bA"  +  h6bB^) (I  +  h6fAn  +  h6fBn) (Qn+1  -  Qn)  = 

y+  X"  y"  (3.10) 

-h  (<5bF+  +  5*f"  +  <5^G+  +  &lG~)n  +  1/3  (Qn  -  Qn-1) 
x  x  y  y 


where  h  =  (2At/3)  and  A+  =  ^q— /  etc.  Equation  (3.10)  can  be 
put  into  its  obvious  algorithm  form  as 


(I  + 

h6bA^  +  h6bB^) AQ*  =  RHS 

(3.11a) 

(I  + 

h6^  +  h6yB”)AQn  =  AQ* 

(3.11b) 

Qn+1 

=  Qn  +  AQn 

(3.11c) 

where  RHS 

represents  the  right  hand  side  of  (3.10). 

The  first 

step  of  the  algorithm  (3.11a)  requires  a  lower  triangular  in¬ 
version  (i.e.  solution)  with  4x4  block  elements.  The  second 
step,  (3.11b),  requires  an  upper  triangular  inversion.  Both 
such  solution  processes  are  simple  compared  to  the  block  tri¬ 
diagonal  inversions  required  with  (2.8)  when  applied  to  only 
inviscid  flow.  The  standard  solution  scheme  (2.8)  is  still 
competitive  with  (3.10),  however,  because  A*  and  B“  are  much 
more  costly  to  form  than  A  and  B. 

Various  other  implicit  algorithms  are  possible  with  flux 
splitting  fc.f.  32],  If  the  thin-layer  viscous  terms  are  in¬ 
cluded  the  following  differencing  has  merit 

[I  +  h ( 6 b A j?  +  5  Bn  +  ?  Mn)][I  +  h5*A?](Qn+1  -  Qn) 

X  *r  y  y  x 

(3.12) 

-  -h(SbF+  +  6*f"  +  6yG  +  IyS)n  +  a (Qn  -  Qn~1)/3 

where  h  ■  —  a  =  0  or  1  for  first  or  second  order  ac¬ 

curacy,  and  6y  and  Ty  are  the  central  difference  operators  de¬ 
fined  previously.  A  solution  algorithm  for  (3.12)  entails 
block  tridiagonal  inversions  carried  on  with  a  forward  sweep 
in  x,  followed  by  a  simple  backsweep  in  the  x-direction. 


The  schemes  given  by  (3.10)  and  (3.12)  have  not  yet  been 
applied  to  as  complex  geometry  situations  as  the  Beam-Warming 
class  of  algorithms  represented  by  (2.8).  The  scheme  (3.10) 
has  been  used  on  stretched  grids  to  compute  inviscid  transonic 
flow  about  a  biconvex  airfoil;  however,  thin  airfoil  boundary 
conditions  were  employed.  A  typical  solution  is  shown  in  Fig. 
12.  This  result  was  computed  without  the  benefit  of  the  tran¬ 
sition  smoothing,  (3.7).  A  viscous  supersonic  wedge  flow  cal¬ 
culation  using  (3.12)  is  indicated  by  Fig.  13.  This  result  is 
an  old  one  that  used  an  earlier  flux  splitting,  namely  (for  F~) 


2X+1 

=  u 

+ 

Ittl 

2X1  - 

u  -  I  u 

+ 

2X3 

=  u 

+ 

1  U 1  +  c 

2X3  = 

u  -  I  u 

2X4 

■  u 

+ 

1  u  1 

2X4  " 

u  -  |  u 

In  this  case,  the  exact  geometry  was  fitted  using  shear  trans¬ 
forms  and  a  very  fine  grid  was  needed  to  resolve  the  viscous 
layer.  A  preliminary  version  of  the  turbulence  model  de¬ 
scribed  in  [9]  was  used  in  the  calculation. 

IV.  REDUCED  SYSTEMS. 

Time-accurate  or  time-like  iterative  methods  are  fre¬ 
quently  used  to  obtain  steady  state  solutions.  If  only  a 
steady  state  solution  is  sought,  however,  one  can  attempt  to 
precondition  and  otherwise  try  to  reduce  the  system  of  partial 
differential  equations  to  obtain  a  more  efficient  solution. 

Not  surprisingly,  certain  reductions  of  the  Euler  equations 
can  begin  to  take  on  features  of  classical  aerodynamic  formu¬ 
lations.  One  such  formulation  [42],  discarded  several  years 
ago  in  favor  of  the  schemes  discussed  earlier,  is  being  re¬ 
vived  because  of  its  excellent  computational  efficiency  for 
steady  rotational  subsonic  flow. 


The  nonconservative  form  of  the  Euler  equations  can  be 
written  in  the  matrix  form 

A3  Q  +  B3  Q  =  0  (4.1) 

a  y 

where 
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If  the  x-axis  is  aligned  with  the  mean  flow  direction,  one  can 
perburb  A  and  B  about  a  reference  state,  o,  with 

u  -  uq,  v  **  0,  p  «  pQ,  p  *  PQ.  (4.2) 

Then  (4.1)  can  be  put  into  the  perturbation  form 

3XQ  +  A^VyQ  *  e  (4.3) 

where 

e  -  (A“1Bo  -  A^B)  3yQ  (4.4a) 


or 

e  -  %Lt(A0-A)3xQ  +  (Bc-B)3yQ]  (4.4b) 

The  e-formulation  (4.4a)  is  obtained  by  multiplying  (4.1)  by 

A"1  and  placing  A~^B  in  the  perturbation  form  K  *  X“*B0  “ 

(27  LB  -  X_1B)  .  To  obtain  the  e-formulation  (4.4b),  A  and  B 
o  o 

are  put  in  perturbation  form  and  the  equations  are  multiplied 

by  27^ .  The  (4.4a)  formulation  has  been  successfully  used  in 
o 

numerical  calculations  (as  described  later)  and  requires  fewer 


derivative  operators  than  the  (4.4b)  formulation.  The  (4.4b) 

formulation  is  not  singular  at  a  sonic  or  stagnation  point, 

but  it  has  not  been  attempted  in  numerical  calculation. 

The  importance  of  (4.3)  is  the  simplicity  of  the  left 

hand  side.  As  indicated  in  (4.2),  v^  is  taken  as  zero  so 

A„  B  has  the  reduced  form 
o  o 


B 


(0  0 

0  0 

0  0 
0  0 


— p  u  ® 
Ho  oyo 

2. 

c 

o  Yo 


-p  u  c  ffl 
o  o  o  yo 


(4.5) 


2  2 

where  =  1/ (c  -  u^  )  and  c  is  the  reference  sound  speed, 
o  o  o  o 

The  eigenvalues  of  a”1Bq  are: 


For  subsonic  mean  flow  the  nonzero  eigenvalues  are  imaginary, 
for  supersonic  mean  flow  they  are  real.  Thus  the  left  hand 
side  part  reflects  elliptic  or  hyperbolic  behavior  as  uq  is 
less  than  or  greater  than  cq. 

As  the  above  matrix  makes  clear,  only  the  third  and 
fourth  equations  of  the  system  (4.3)  are  strongly  coupled. 
Written  out  we  obtain 


(4.6a) 

(4.6b) 

(4.6c) 


(4 . 6d) 


Thus  if  are  known  from  some  previous  estimate,  we  can  solve 
equations  (4.6c)  and  (4.6d)  for  v  and  p.  Once  v  is  obtained, 
p  and  u  are  obtained  from  (4.6a)  and  (4.6b)  via  simple  inte¬ 
gration.  Better  estimates  of  e.  can  now  be  formed  and  the 
process  is  repeated  until  iterative  solution  of  (4.1)  is  ob¬ 
tained. 

The  equations  (4.6c)  and  (4.6d)  form  an  elliptic  system 
when  uq  <  cq.  This  is  readily  apparant  as  the  eigenvalues  of 


-p  u  c  9 
o  o  o  o 


(4.7) 


are  imaginary  as  indicated  previously.  The  equations,  which 
if  linearized  by  setting  =  0  could  be  transformed  into 
Cauchy  Riemann  equations,  can  be  differentiated  to  form  a 
Poisson  equation  with  either  v  or  p  as  the  dependent  variable. 
For  example 

tl-<o0/-=o|2|3«v  -  Syyv  -  Hl-(u0/c0)2!  (Vj-'VVSs* 


Once  v  is  obtained,  p  is  found  by  integrating  (4.6d).  If 
uQ  >  cQ,  (4.8)  represents  a  wave  equation  and  the  eigenvalues 
of  (4.7)  are  real. 

The  (4.4a)  formulation  using  (4.8)  to  update  v  has  been 
successfully  applied  in  two-dimensions  [42]  to  compute  rota¬ 
tional  subsonic  flow.  Typical  results  for  a  lifting  biconvex 
airfoil  are  shown  in  Fig.  14  for  the  incoming  shear  flow  (jet 
or  defect)  profiles  as  depicted  in  Fig.  15.  Computational 
times  for  a  fully  converged  solution  on  a  59x90  grid  averaged 
18  seconds  per  case  on  a  Control  Data  7600  computer. 

The  algorithm  using  (4.8)  is  robust.  For  inexplicable 
reasons,  elimination  of  v  from  (4.6c)  and  (4.6d)  and  solution 
of  a  Poisson  equation  for  pressure,  p,  has  always  been  a  dis¬ 
appointment.  The  failure  is  believed  to  be  keyed  to  the 
boundary  condition  treatment. 


The  above  ideas  readily  extend  to  three-dimensions  al¬ 
though  no  calculations  have  yet  been  undertaken.  The  three- 
dimensional  perturbation  form  of  the  equations  is  given  by: 


(4.9) 


In  this  case,  v,  w,  and  p  are  strongly  coupled  through 


-1, 


3„v  +  (pQuo)  3„p  =  e 


3xw  +  (pouo)"l3zp  "  e 


3xp  “(pouoco  *oJ  (3yv  +  3zw)  *  e5 


(4.10a) 

(4.10b) 

(4.10c) 


Once  v,  w,  p  are  obtained,  p  and  u  are  found  by  integration  of 


3xp  "  £1  +  (poVo}  (3yv  +  3zw) 


3xu  *  e2  -  (co  V  (3yv  +  3zw) 


(4 . lOd) 

(4 . lOe) 


Equations  (4.10a)  to  (4.10c)  can  be  differentiated  so 
that  v  and  w  are  eliminated  to  form  a  Poisson  equation  in 
pressure,  that  is 


tl-<Vco>  ]3xxp  +  3yyp  +  3zzp  "  * 


(4.11) 


Alternately,  as  this  was  unsuccessful  in  two-dimensions,  pres¬ 
sure  can  be  eliminated  and  vector  potential  like  equations  can 
be  formed  for  v  and  w.  In  particular 


tl-<uo/co)‘]3xxv  ♦  *yyV  +  3yzw  -  f3 


(4.12a) 


[l-(u0/c0)  13xxw  +  3zzw  +  3yzv  -  f. 


(4.12b) 


Here  again,  if  (4.12a)  and  (4.12b)  are  solved  for  v  and  w, 
then  p,  p,  and  u  can  be  found  by  simple  integrations  of 
(4.10c)  to  (4.10e) . 

Although  not  discussed  previously,  one  can  draw  some  in¬ 
teresting  analogies  between  (4.10a)  -  (4.10c)  to  the  incom¬ 
pressible  irrotational  relations 


3xu  +  3yv  +  3zw  =*  0  (4.13a) 

3  w  -  3  v  «  0  (4.13b) 

y  z 

3  u  -  3  w  *  0  (4.13c) 

Z  X 

3  v  -  3  u  *  0  (4.13d) 

x  y 


The  potential 

4>x  -  u,  -  v,  -  w  (4.14) 

satisfies  (4.13b)  to  (4.13d)  and  from  (4.13a)  gives  the 
Laplacian 


$  +  $  +  <ji  ■  0 

xx  yy  zz 


(4.15) 


The  particular  vector  potentials  x  defined  as  [43] 


u  *  ij>  +  X- *  v  *  ,  w  ■  -x„ 

y  z  x  x 


(4.16) 


satisfies  (4.13a)  and  from  (4.13c)  and  (4.13d)  give 


ill  +  il>  +  x  *  0 
rxx  ryy  Ayz 


(4.17a) 


X  +  X  + 
Axx  Azz 


yz 


o 


(4.17b) 


The  analogies  between  (4.11)  and  (4.15)  and  between 
(4.12)  and  (4.17)  are  completely  transparent  and  extend  even 
to  the  boundary  condition  treatment.  Curiously,  (4.15)  is 
elliptic,  but  (4.17)  is  not  strictly  elliptic  according  to  the 
definition  of  [44].  Moreover,  the  system  of  equations  com¬ 
prised  of  (4.13a),  (4.13c),  and  (4.13d)  is  not  elliptic.  This 
latter  system  presents  no  difficulty  to  numerical  solution, 
however,  and  conventional  numerical  algorithms  such  as  succes¬ 
sive  overrelaxation  have  been  successfully  applied  to  the  sol¬ 
ution  of  (4.17).  It  is  conjectured,  therefore,  that  one  could 
devise  rapid  solution  procedures  for  (4.12)  and  (4.10c) 
through  (4.10e)  since  a  similar  approach  worked  in  two-dimen¬ 
sions. 

The  perturbation  schemes  described  here  have  only  been 
used  in  numerical  calculation  of  pure  subsonic  two-dimensional 
flow.  Extensions  to  pure  supersonic  flow  appear  to  be 
straight  forward  but  such  is  not  the  case  for  transonic  flow. 
The  isentropic  primitive  variable  approach  of  Martin  [45]  and 
the  stream  function  sonic  line  treatment  of  Hafez  [46]  may  of¬ 
fer  guidance,  however.  Note  there  is  no  problem  bringing  the 
essentials  of  the  perturbation  formulation  into  conservative 
form.  For  example 

3  ,5  *  ^Vy5  *  Soll3*lV-BolF)+3yl5o5-“olGI  1  l4-18) 


where  the  matrix  N  relates  Q  to  Q  via 

N  - 

3  Q 

and  Q,  F,  G  are  the  Cartesian  form  of  the  conservative  fluxes 
of  (2.1).  The  left  hand  side  of  (4.18)  is  identical  to  (4.3). 

Rather  than  attempt  to  extend  the  method  of  this  section 
into  the  realm  of  transonic  flow,  it  may  be  better  to  use  this 
procedure  as  the  outer  part  of  a  chimera  (i.e.,  multiple  con¬ 
stituent)  algorithm.  For  example,  the  semi-flux-split  scheme 
(3.12)  can  be  used  to  solve  supersonic  and  subsonic  flow,  and 
it  is  very  efficient  for  supersonic  and  high  subsonic  regions. 


The  scheme  of  this  section  is  very  efficient  for  subsonic 
flow.  A  combined  chimera  method  for  transonic  flow  simulation 
about  an  airfoil  might  entail  use  of  the  schemes  (4.3)  and 
(3.12).  The  outer  subsonic  flow  would  be  solved  with  (4.3), 
the  inner  high  subsonic  and  supersonic  flow  with  (3.12),  and 
both  solution  regions  could  be  slightly  overlapped  [42,47]  for 
iterative  efficiency. 

V.  CONCLUDING  REMARKS. 

The  successful  simulation  of  the  flow  about  really  com¬ 
plex  geometries  will  require  further  advances  than  those  dis¬ 
cussed  here.  Nevertheless,  the  combination  of  implicit  fi¬ 
nite  difference  procedures,  generalized  coordinates,  and  nu¬ 
merical  grid  generation  techniques  is  proving  to  be  effective. 
The  overall  methodology  can  be  built  on,  and  improvements  in 
computational  efficiency  can  be  expected.  In  this  overall 
development  it  would  seem  that  the  next  crucial  step  is  to 
develop  efficient  chimera  schemes  -  that  is,  modular ly  coded 
numerical  schemes  that  blend  or  overlay  more  than  one  grid 
system  and  more  than  one  type  of  governing  equation  or  numeri¬ 
cal  algorithm. 
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Sketch  1.  Showing  patched  grid  system  with  all  body  surfaces 
mapped  to  grid  outer  boundaries.  This  approach  requires  in¬ 
terpolation  along  boundaries. 


Sketch  2.  Showing  overset  grid  system.  In  this  approach  var 
ious  interior  points  must  be  turned-off  and  interpolation  is 
required  at  grid  interfaces,  however,  grids  are  easy  to  gener 
ate  and  have  minimum  distortion. 


Fig.  1.  Cascade  grid  generated  using  elliptic  partial  differ¬ 
ential  equations. 


Fig.  3.  Schematic  defining  gi 
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Fig.  4.  Computed  variation  of  a 
NACA  65  213  airfoil. 
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Pig.  7.  Mach  contour  lines  during  an  aileron  buzz  cycle, 
0.82,  a=-l  deg,  R  =20.7xl06.  Here  9=(t/t  )*  360  deg  where  t 

6  t'  r 

is  the  time  to  complete  one  cycle. 


On  Application  of  Body  Conforming  Curvilinear  Grids 
for  Finite  Difference  Solution  of  External  Flow 
in  Numerical  Grid  Generation 


from 


Numerical  Grid  Generation 
J.  F.  Thompson,  Editor 
North-Holland,  1982 


ON  APPLICATION  OF  BODY  CONFORMING  CURVILINEAR  GRIDS 


FOR  FINITE  DIFFERENCE  SOLUTION  OF  EXTERNAL  FLOW 
Joseph  L.  Steger 

Department  of  Aeronautics  and  Astronautics 
Stanford  University,  Stanford,  CA  04305 

I.  INTRODUCTION 

Finite  difference  practitioners  frequently  make  use  of  arbitrary  coordinate  transforms 
and  introduce  body  conforming  curvilinear  grid  systems.  The  coordinate  transforms  may 
either  be  built  in  globally  in  mappings  from  physical  space  to  computational  space,  or 
they  may  be  built  in  locally  in  the  finite  volume  sense.  The  advantages  of  using  body 
conforming  curvilinear  grids  in  finite  difference  flow  field  simulation  include  the  following: 
Body  conforming  gTids  simplify  the  application  of  boundary  conditions  insofar  that  grid 
lines  will  coincide  with  the  body  boundary.  Curvilinear  grids  may  be  clustered  to  Bow 
field  action  regions  to  improve  solution  accuracy.  Body  conforming  grids  may  allow 
simplification  of  th'  .  .rning  equations.  Such  grids  can  also  help  maintain  a  well-ordered 
system  of  algebraic  equations  suitable  for  vector-computer  processing  or  approximate- 
factorization-implicit  techniques. 

The  task  of  generating  suitable  body  conforming  curvilinear  grids  is  not  an  easy  one. 
The  grids  should  be  generated  in  an  automatic  manner  requiring  minimum  user  input. 
Yet  the  user  will  wish  to  maintain  considerable  control  of  where  points  will  be  distributed 
along  the  boundary  surface  and  how  they  are  clustered  in  the  interior  field.  Moreover  the 
grid  must  be  tailored  in  some  degree  to  the  numerical  algorithm  because  some  numerical 
algorithms  are  more  sensitive  than  others  to  grid  smoothness,  skewness,  and  stretching. 

Although  use  of  body  conforming  curvilinear  grids  can  offer  the  advantages  cited 
previously  their  careless  application  can  lead  to  difficulties.  This  is  particularly  true  when 
the  governing  equations  arc  differenced  in  conservative  (i.e.,  divergence)  form  and  transform 
metric  terms  are  brought  inside  the  difference  operators.  Then  as  noted  above,  some 
numerical  algorithms  are  far  more  mesh-sensitive  than  others,  and  numerical  accuracy 
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j  and  computational  efficiency  can  be  affected  by  bow  rapidly  a  grid  changes  or  bow  far 
3  away  it  is  from  orthogonality. 

_  The  subject  of  this  paper  is  not  the  generation  of  body  conforming  curvilinear  grids; 
~  rather  it  is  the  use  of  such  grids  in  finite  difference  applications.  Jn  Section  D  of  this  paper 
~  the  difficulties  of  solving  the  transformed  equations  in  conservative  form  are  discussed. 

-  In  Section  III  various  experiences  are  cited  to  suggest  that  considerable  computational 

-  efficiency  can  yet  be  gleaned  by  further  improvements  of  the  grid.  Concluding  remarks 
__  follow  in  Section  IV. 

1 

]  H.  CONSERVATIVE  DIFFERENCING  OF  TRANSFORMED  EQUATIONS 

1 

—  a)  Background 

In  aerodynamics  applications  we  frequently  try  to  difference  the  governing  equations 
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in  conservative  or  divergence  form.  Conservative  form  differencing  is  preferred  because 
it  best  maintains  the  correct  weak  solution  of  the  governing  equations.  Thus  if  a  shock 
wave  is  captured  by  simply  solving  the  difference  equations  (as  opposed  to  fitting  a  shock 
wave  discontinuity  into  the  difference  equations),  then  the  speed,  location  and  jump  of 
the  shock  can  only  be  correct  if  the  partial  differential  equations  are  in  conservative 
‘  form.  The  difference  equations  must  also  satisfy  the  divergence  relation,  at  least  in  the 

—  vicinity  of  the  shock.  The  conservative  form  of  the  equations  may  also  be  desirable  for  its 

-  numerical  properties.  For  example,  nonlinear  equations  in  conservative  form  can  be  more 
J  cleanly  linearized  about  a  previous  state  than  those  in  nonconservative  form.  This  can 

~be  advantageous  iu  implicit  marching  procedures  in  order  to  avoid  iterative  solutions  of 
~  nonlinear  equations  with  each  marching  step. 

Let  the  conservat ion-law- form  of  the  equations  be  represented  in  Cartesian  coordinates 


+  dyG  —  0 


lines  to  ~7 
BOTTOM  - 

e 


s  where  for  simplicity  only  two  dimensions  are  considered.  This  strict  divergence  form  of 
? 

j  the  equations  can  be  maintained  for  new  independent  variables 
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MICROCOPY  RESOLUTION  TEST  CHART 

NATIONAL  3URtAU  OF  STANDARDS  1963  A 


£  =  £(*.  it,  0 

n  =  *(*,  v,  0 

t  *=  t 


(2) 


« M-  U|) 


9r^  +  =  0 


where 


£’“^-l(£.«  +  £.£,  +  £,C) 

G  *=  +  ij*F  +  ntC) 


end  where  J  is  the  transform  icobian 


J  —  6f»  -  £»*» 


(3) 

(4o) 

(46) 

(*c) 

(5) 


For  a  thermally  perfect  gas  Q,  F,  and  G  may  represent  the  Cartesian  inriscid  and 
viscous  dux  quantities  for  conservation  of  mass,  momentum  and  energy.  For  example,  for 
inviscid  flow 


Q  = 


f) 

pu 

pv 

,  r  = 

/  p»  > 

P**+P 
pu  V 

,  c  = 

/  i 

pav 

pv*  +P 

v  J 

+  P)> 

U«+p)> 

(6) 


where  p  is  fluid  density,  u,  v  and  to  are  Cartesian  velocities  components,  p  is  pressure  and 
e  is  given  by 


«-(h  rlp+\f(«*+v*) 


(7) 


Alternately  in  the  case  of  compressible  potential  flow 


Q=*p,  F=*p4„  G=*p<j>t  (8) 

and  p  =  p{<t>)  is  determined  by  the  Bernoulli  relation. 

Although  the  transformed  governing  equations  (3)  are  more  complex  than  (1),  apparent 
simplicity  is  returned  to  the  inviscid  flow  equations  with  introduction  of  the  unsealed 

3 


All  such  right-hand-side  combinations  of  metric  terms  are  found  to  be  zero  because  of  the 
relations 


v, 

ijj  -  -x, 

ft  =*  ~Zr(m  -  Jf,ft 


if,//  =  -y( 

Ifi  *=>  -  »r»?» 


for  example, 


F{V(n  -  »,<)  »  0 


(IS) 


(16) 


and  so  on. 

The  point  of  ail  this  is  that  the  metric  quantities  have  been  worked  inside  the  differential 
operators.  This  is  possible  because  combinations  of  metric  terms  such  as 


40)+8'(~O  =  ° 


(17) 


are  found.  Note  also  that  unlike  equation  (1),  equation  (3)  has  been  scaled  by  /“*. 

b)  Metric  Differencing 

The  fact  that  the  transformed  governing  equations  now  have  the  metrics  brought 
inside  the  difference  operations  can  lead  to  numerical  errors.  This  b  because  the  metric 
variation  is  now  being  differenced  along  with  the  Dow  field  quantity.  In  typical  external 
aerodynamics  applications,  the  flow  quantities  far  from  the  body  are  essentially  constant 
or  uniform.  Difference  terms  should  therefore  be  zero,  and  this  will  be  true  if  equation  (1) 
la  differenced  on  a  uniform  mesh.  For  a  nonuniform  mesh  the  transformed  equations  will 
not  yield  zero  in  regions  of  constant  flow,  however,  unless  proper  differences  of  the  metric 
identities  are  zero,  that  is  from  (M) 


*t(ft/^  +  <e(ff.//)-0 


(««) 


M«r//)  +  «e(*AO-0  (1«4) 

+  +  (18c) 

where  ft,  ft  and  ft  are  the  difference  operators  used  in  the  solution  algorithm  for  (3).  If, 
for  example,  the  metrics  ft//,  etc.,  could  be  exactly  evaluated  (as  they  can  be  in,  say, 
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•  cylindrical  coordinate),  then  equations  (18)  will  not  be  exactly  zero  but  will  be  zero 
to  within  the  order  of  accuracy  of  the  operators  S(,  etc.  For  a  rapid  variation  of  the 
metrics  and  for  large  grid  spacing — a  phenomenon  frequently  occurring  in  aerodynamic 
applications  in  the  far  field — this  error  can  be  very  appreciable.  It  can,  in  fact,  add  a 
error  source  to  the  equations  that  can  overwhelm  the  solution  accuracy.  However,  if  the 
metrics  themselves  are  differenced  such  that  equations  (18)  are  exactly  zero,  then  this  error 
is  controlled.  For  example,  in  two  dimensional  steady  flow  relations  (18a)  and  (18b)  become 
■sing  (15) 

*<(*,)  -«,(»)  “0  (19«) 

-dt(r„) +  *„(*<)  =  0  (104) 

If  and  are  differenced  as  Sny  and  where  Sn  and  5;  are  the  same  difference  operators 
used  in  the  solution  algorithm  for  (3),  then  the  metric  identities  (19)  exactly  difference  to 
zero.  The  importance  of  satisfying  these  relations  was  pointed  out  in  [2]  in  which  three 
point  central  spatial  differences  were  used  in  the  solution  algorithm  [3]  as  well  as  for  the 
metric  quantities. 

In  three  dimensions  it  becomes  more  difficult  to  exactly  difference  the  metric  identity 
relations.  For  steady  or  simple  unsteady  grid  motion,  Pulliam  and  Steger  (4)  introduced 
an  averaging  process  for  the  steady  terms  that  works  for  any  difference  operator  that  can 
be  differenced  in  parts  as 

S[nv)  =  Oi»K«u)  +  (/i«K*»)  (W) 

where  p  is  an  averaging  operator.  An  example  of  (20)  is  given  by 

Vu»»*(fiv)(Vu)  +  (i»u)Vt>  (21) 

where  ?•■«/-  uy_t,  pu  *  ?l  ,  etc.  In  extended  work  Thomas  and  Lombard  [5] 
correctly  treated  the  unsteady  metric  variations  and  cleverly  simplified  the  calculation  of 
the  special  metric  terms.  They  also  coined  the  term  "metric  conservation  law”  to  describe 
the  fact  that  the  metric  relations  (18)  must  be  differenced  to  be  zero. 

As  one  introduces  higher  order  central  or  one-sided  spatial  difference  operators,  or 
uses  predictor-corrector  schemes,  it  becomes  more  and  more  difficult  to  correctly  satisfy 
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the  metric  relations  (18).  This  ultimately  lead  (4|  to  an  approximate  caoceilatioD  method 
that  relies  on  solving  the  differenced  equations  in  a  simple  perturbation  form: 

MS-O +«<(*•-£») -m,(6-6»)=*o  (22) 

In  the  far  field  etc.,  and  any  consistent  difference  scheme  is  satisfied  regardless  of 

•  how  rapidly  the  metrics  vary.  When  F  appreciably  varies  from  fee,  etc.,  it  is  assumed  that 
the  grid  is  sufficiently  smooth  and  refined  so  that  the  metric  error  is  no  greater  than  the 
error  of  differencing  the  flux  terms.  In  external  flow  applications  this  process  has  worked 
quite  well,  including  successful  use  with  the  potential  flow  equation  (cf.  (6|). 

A  direct  means  of  cancelling  the  metric  errors  is  the  straightforward  one  of  subtracting 
the  error,  for  example,  for  a  stationary  grid 

StQ  +  Sft"  +  *■  F(S(  jr,  —  Sny()  +  G{—S(Xn  +  Snt()  (23) 

'  This  puts  the  equation  into  a  weak  conservation  law  form  that  in  principle  does  not  degrade 
the  shock  capturing  properties  of  the  scheme.  It  could,  however,  contribute  to  a  mild 
source-term  weak  instability  that  would  be  alleviated  somewhat  by  spatial  averaging  of 
the  right-hand  side  F  and  G  flux  terms.  These  right-hand-side  fiux  terms  can  also  be 
approximately  evaluated  at  oo  in  somewhat  duplication  of  (22)  above. 

The  whole  problem  of  differencing  the  metrics  has  been  avoided  in  [7-IOj.  In  this 
approach  the  Cartesian  equations  are  expanded  by  chain  rule  and  then  simply  left  that 
way.  That  is 

9rQ  +  (td(Q  +  ltd, ,Q  +  £ *dfF  +  t\xdnF  +  «=  0  (24) 

and  is  called  the  quasi-linear  form  by  Sbamortb  and  Gibeling  [10]  or  chain-rule  conservation 
form  by  Hindman  |9].  Although  the  Jacobian  is  never  divided  through,  this  form  is 
somewhat  similar  to  the  weak  conservation  law  form  (23)  just  discussed,  particularly  so 
with  averaging  of  F  and  G.  It  should  also  properly  capture  the  correct  jump  relations. 
The  chain-rule  conservation  form  may  well  be  a  good  compromise  to  differencing  the 
transformed  equations  in  conservative  form.  For  certain  algorithms  (e.g.,  Beam-Warming 
(3|)  it  appears  to  require  more  work  than  using  the  strong  conservation  law  form  with  free 
stream  correction. 
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e)  Perturbation  Form  Digression  * 

The  above  idea  of  subtracting  out  the  free  stream  metric  variation,  equation  (22), 
discussed  previously  is  a  special  ease  of  perturbing  the  solution  about  a  known  function 
which  in  some  sense  is  also  a  nearby  or  approximate  solution.  Let  Q  =  Qo  +  O'  where 
Qt  is  the  nearby  or  approximate  solution  and  let  O'  be  the  perturbation.  The  terms  of 
equation  (3)  can  be  rewritten  as, 


d,Q  —  d,Q0  +  d,Q' 
d(t'lQ)  -  9(F(Q0)  +  d((t{Q)  -  /•(<?•)) 

etc. 

Then  assuming  functions  of  Qo  are  sufficiently  simple  to  be  very  accurately  (or  exactly 
differentiated)  with  operators  S,  the  differencing  of  equation  (3)  can  be  represented  as 


6rQ’  +»( 


[*(<?)  -  F(Q0)j  +  *,[6(Q) "  6(<?o)] 


(25) 


where  ST,  S(,  and  5,  represent  the  algorithm  difference  operators  and  S„  S(,  represent 
the  very  accurate  differencing.  In  the  case  Qo  **  Q*>  the  right-band  side  is  analytically 
xero  and  equation  (25)  is  identical  to  equation  (22). 

Such  a  perturbation  form  of  the  differenced  equation  has  been  proposed  in  internal 
spin-up  problems  to  remove  the  ax  {symmetric  variation  from  the  dependent  Cartesian 
velocity  variable  (11,12).  It  might  also  be  used  in  problems  in  which  certain  fine  details 
might  be  otherwise  lost  in  a  coarse  grid.  For  example,  a  nonuniform  incoming  flow  profile 
can  be  represented  in  Qo  that  would  otherwise  be  lost  in  a  far  field  coarse  grid.  Assuming 
in  this  case  that  Qo  satisfies  the  Euler  equations,  the  right-hand  side  of  (25)  a  identically 
xero.  Calculations  using  this  particular  technique  for  incoming  inviscid  shear  Sows  have 
been  tested  by  Buning  and  Steger  [13],  Although  not  yet  tried,  Qo  might  be  chosen  as 
an  approximate  solution,  in  which  case  the  right  hand  side  of  (25)  would  not  be  xero.  In 
regions  in  which  Q-*Qo,  one  could  hope  to  use  a  much  coarser  grid  without  losing  solution 
accuracy. 
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HI.  CALCULATIONS  ON  CURVILINEAR  GRIDS 

Finite  difference  and  related  finite  volume  calculations  using  body  fitted  curvilinear 
grid  systems  have  been  carried  out  for  some  time.  Solution  variables  have  included 
velocity  potential,  stream  function,  and  the  primitive  variables.  Computed  results  include 
incompressible  and  compressible  Bow  around  airfoils,  projectiles,  cascades,  inlets,  wings, 
wing-body  combinations,  etc.  In  some  cases  the  body  deforms  with  time  (e.g.,  airfoil 
with  moving  aileron)  and  occasionally  solutions  for  multiple  non-connected  bodies  appear 
(e.g.,  airfoil  with  detached  flap).  No  attempt  will  be  made  to  review  this  work — the 
interested  reader  will  find  much  of  the  materia)  in  the  AIAA  Journal,  the  Proceedings  of 
the  International  Conference  on  Numerical  Methods  in  Fluid  Mechanics,  Computers  and 
Fluids,  and  the  Journal  of  Computational  Physics.  What  b  apparent  from  thb  literature  b 
that  while  we  are  becoming  more  adept  at  solving  the  flow  about  complex  configurations, 
considerable  computational  efficiencies  are  yet  to  be  obtained. 

In  order  to  illustrate  points  to  be  discussed  later,  the  results  of  a  finite  difference 
simulation,  due  to  Nicolet  et  at.  [ldj  for  flow  about  an  X-24C  configuration  b  reproduced 
in  Figs.  1  to  3.  A  head-on  view  of  the  X-24C  b  indicated  in  Fig.  I,  while  Fig.  2  shows 
typical  views  of  the  grid  fit  between  the  body  surface  and  an  analytically  fit  outer  bow 
shock.  The  grid  in  thb  case  b  generated  with  a  hyperbolic  partial  differential  equation  grid 
generation  scheme  [1S|.  The  overall  three  dimensioaal  grid  b  formed  by  generating  two 
dimensional  grids  at  each  station  along  the  body  as  the  solution  progresses  by  marching 
the  steady  parabolised  Navier-Stokes  (PNS)  equations.  The  hyperbolic  grid  generation 
scheme  b  fast  enough  to  be  used  within  the  fiow  Add  marching  scheme.  Moreover,  each 
two  dimensional  grid  b  itself  generated  using  the  same  kind  of  numerical  algorithm  used 
for  solving  the  PNS  equations— a  sort  of  conservation  of  numerical  algorithm  knowledge. 
Computed  surface  pressure  and  heat  transfer  at  a  station  just  prior  to  the  beginning  of 
the  wing  are  compared  to  experimental  data  in  Fig.  3. 

In  carrying  out  the  preceding  calculation  or  any  similar  calculation  on  a  generalised 
grid  it  b  found  that  the  solution  accuracy  depends  'on  the  grid.  Thb  b  not  surprbing 
because  unless  one  has  a  very  fine  mesh  throughout  the  field,  an  accurate  solution  will 
require  that  grid  points  be  clustered  to  the  flow  field  action  regions— the  change  of  gradient 
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regions.  The  grid  in  Fig.  2,  for  example,  is  exceedingly  fine  near  the  body  surface  in  order 
to  resolve  viscous  gradients  along  the  wall.  (The  less  than  perfect  agreement  with  the 
experimental  heating  rate  shown  in  Fig.  3  is  apparently  not  due  to  inadequate  resolution 
in  this  direction).  The  outer  grid  line  also  coincides  with  the  bow  shock,  and  points  are 
clustered  along  the  body,  for  example  to  resolve  the  cross  flow  expansion  around  the  chime 
(i.e.,  lower  right  corner  in  Fig.  2)  when  the  vehicle  b  at  angle  of  attack  (here  at  ft*). 
Otherwise  no  other  attempt  was  made  in  this  calculation  to  adapt  the  grid  to  computed 
flow  field  gradients. 

Besides  proper  grid  clustering,  the  smoothness  of  the  grid,  the  skewness  of  the  grid, 
and  sometimes  the  aspect  ratio  of  the  grid  can  affect  the  accuracy  of  a  numerical  solution 
or  the  efficiency  with  which  it  is  obtained.  The  grids  shown  in  Fig.  2  are  nearly  orthogonal 
close  to  the  body  surface  and  they  have  a  smooth,  gradual  variation.  These  grids  would 
be  judged  felicitously.  However,  the  quality  of  a  grid  seems  to  be  hard  to  quantify  be¬ 
cause  various  numerical  algorithms  appear  to  behave  differently  to  the  properties  of  grid 
smoothness,  skewness,  and  stretching.  Numerical  algorithms  that  use  a  very  compact 
stencil  of  points  to  evaluate  fluxes  and  metrics,  for  example,  generally  seem  to  be  less  sen* 
sitive  to  grid  spacings  that  change  rapidly  or  even  discontinuously.  Thus  computors  using 
the  MacCormack  finite  volume  method  for  Navier-Stokes  equations  sometimes  change  the 
grid  spacing  by  a  factor  or  2  or  4  in  a  given  region.  Such  a  change  b  not  allowed,  for 
example,  when  using  high  order  central  spaeial  differencing  operators. 

Some  numerical  algorithms  appear  very  sensitive  to  mesh  cell  aspect  ratio,  Le.,  the 
ratio  of  Ax  to  Ay  or  (*$  +  yf)1^*  to  (x*  +  Thus  Jameson  [16J  in  developing 

a  multtxrid  relaxation  algorithm  for  the  transonic  potential  equation  abandoned  SLOR 
iterative  methods.  In  its  place  he  used  an  aIternating-direction*impiicit  scheme  as  the 
multigrid  iterative  solver  because  it  b  less  sensitive  to  cell  aspect  ratio.  The  very  efficient 
approximate-factorization-implicit  relaxation  scheme  of  Holst  [17|  appears  to  degrade  if 
uniform  fine  grid  spacing  b  used  along  the  body,  prompting  Hobt  to  generate  hb  grids 
with  tbb  constraint  in  mind.  For  hb  approximate-factorization  scheme  the  grid  shown  in 
Fig.  4  b  much  preferred  to  that  shown  in  Fig.  S.  A  user  of  a  standard  alternating-direction- 
implicit  relaxation  scheme,  however,  may  very  well  opt  instead  for  the  grid  of  Fig.  S  over 
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that  of  Pig.  4  simply  because  of  its  finer  grid  spacing  near  the  body  and  presumably  greater 
accuracy. 

Avoiding  certain  undesirable  grid  properties  such  as  skewness  can  lead  to  more  com¬ 
plex  computer  programs  and  perhaps  other  difficulties.  The  cascade  C-grid  shown  in  Fig. 
ft  for  example  is  too  highly  skewed.  While  the  implicit  Beam- Warming  algorithm  for  the 
Euler  or  Navier-Stokes  equations  functions  on  such  a  grid,  it  runs  far  from  optimum.  An 
alternative  grid  to  that  shown  in  Fig.  8  might  use  an  overlapped  or  patched  grid  system. 
For  example,  the  overlapped  grid  system  shown  schematically  in  Fig.  7  b  suggested  because 
each  grid  b  easy  to  generate  and  has  minimum  dbtortion.  However,  such  a  grid  system 
requires  extensive  modification  of  exbting  numerical  algorithms  and  computer  programs. 
Thb  b  because  certain  grid  points  will  have  to  be  flagged  ofT,  and  grid  interfaces  will  have  to 
be  joined  without  degrading  numerical  stability.  Nevertheless,  overlapped  or  patched  grid 
'  systems  will  ultimately  be  needed  as  body  boundary  configurations  become  more  complex, 
for  example,  in  computing  flow  about  a  wing  with  engine  nacelles. 

Finally,  it  should  be  remarked  that  the  effect  of  a  poorly  spaced  grid  will  sometimes 
not  be  observed  until  the  data  b  displayed  or  utilized  in  a  different  way.  The  unpublbhed 
result  due  to  Seidel  [18]  that  b  shown  in  Fig.  8  is  an  example.  The  plots  of  generalized 
pitching  moment  versus  reduced  frequency  show  an  essentially  exact  solution  (dashed  line) 
and  a  low  frequency  transonic  small  dbturSance  finite  difference  solution  with  the  nonlinear 
terms  removed.  The  flow  b  about  a  flat  plate  subjected  to  an  angle  of  attack  pulse.  The 
amall  oscillations  shown  in  the  finite  difference  result  are  a  significant  error  in  a  flutter 
calculation.  They  were  traced  back  to  a  discontinuous  change  of  grid  stretching  more 
than  a  chord  away  from  the  airfoil,  and  were  eliminated  by  using  a  smoothly  stretched 
grid  throughout. 

IV.  CONCLUDING  REMARKS 

Finite  difference  methods  coupled  with  body  conforming  curvilinear  grid  systems  are 
being  used  to  solve  a  variety  of  complex  flow  fields.  Current  numerical  algorithms  anti 
grids  are  tuned  to  flow  field  applications  that  can  be  computed  in  reasonable  times  on 


present  day  machines.  These  numerical  algorithms  rely  on  sparse- equation  time-accurate 
or  iterative-solution  methods  that  function  best  on  well  ordered  grids. 

Curvilinear  body  conforming  grids  have  made  modern  finite  difference  schemes  into 
practical  engineering  tools.  They  simplify  the  application  of  boundary  conditions  and  allow 
flow  field  gradients  to  be  resolved  in  an  orderly  manner.  However,  one  must  be  careful  in 
differencing  transformed  equations,  especially  when  the  equations  are  in  conservative  form 
and  transform  metrics  are  brought  inside  the  difference  operators.  Finite  difference  algo¬ 
rithms  are  also  sensitive  to  grid  smoothness,  skewness  and  stretching  with  some  algorithms 
being  much  more  adversely  affected  than  others. 

As  finite  difference  methods  are  extended  to  more  complex  geometries,  it  becomes 
obvious  that  more  than  one  grid  system  will  have  to  be  used.  Exactly  how  multiple 
grids  should  best  be  joined,  patched,  or  overset  together  remains  a  research  topic,  but  a 
.number  of  approaches  will  likely  give  satisfactory  results.  The  simultaneous  development  of 
multiple  grid  systems  and  finite  difference  schemes  suitable  for  multiple  grids  is  underway 
and  will  be  a  major  pacing  item  in  computational  fluid  dynamics. 
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Fig.  7  Schematic  indicating  how  overlapped  grids  could  be  used  for  cascade  Bow  applica¬ 
tions  to  remove  grid  distortion. 
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Fig.  8  Generalised  pitching  moment  versus  reduced  frequence  showing  oscillation*  due  to 
a  poorly  stretched  grid. 


V.  FUTURE  WORK 


Two  m*w  flow  simulation  codes  for  projectiles  are  now  underway  at  BRL.  One 
of  t  hese,  which  is  being  implemented  by  Nietubicz,  extends  our  segmented  projectile 
code  for  the  calculation  of  ringed  or  tubular  projectiles.  This  advanced  projectile 
flies  a  more  level  trajectory  than  conventional  rounds.  Figures  1  and  2  show  a 
typical  configuration  and  a  preliminary  result  obtained  by  Nietubicz. 

The  other  new  projectile  code  is  still  in  the  planning  stage.  A  three  dimensional 
transonic  projectile  code  with  base  will  use  the  same  segmentation  process  used 
with  the  two  dimensional  code.  In  order  to  minimize  the  number  of  grid  points 
needed  to  resolve  a  three  dimensional  field,  we  plan  to  use  a  spectral  method  in  the 
circumferential  direction.  Although  our  solution  method  is  implicit,  the  spectral 
method  only  has  to  be  implemented  explicitly.  (This  has  been  shown  by  K.  C. 
Reddy  and  has  been  independently  verified  by  my  student,  Mr.  T.  Barth.)  This 
code  should  be  operational  by  December  1984.  My  involvement  will  continue  as  an 
employee  of  NASA  Ames. 
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